REMERCIEMENTS

Remis de mes émotions quelques jours après l'obtention du grade de docteur de l' École Nationale Supérieure des Arts et Métiers, voici venu le temps de remercier les différentes personnes que j'ai pu rencontrer au cours de cette thèse, mais aussi celles qui m'ont accompagné tout au long de cette dernière et bien évidemment celles sans qui je n'aurais jamais pu mener à bien ce projet, qui fut un travail de longue haleine.

Je tiens à remercier dans un premier temps M. Laurent PIZZAGALLI, Directeur de recherche CNRS à l'Institut Pprime de Poitiers, d'avoir présidé ce jury de thèse ainsi que pour les questions et discussions qui ont suivi la synthèse de mes travaux.

Je remercie sincèrement M. Patrick CORDIER, Professeur à l'Université de Lille 1, d'avoir accepté la lourde tâche de rapporteur. Je le remercie également pour ses commentaires et remarques qui ont été très instructifs et dont je me servirai sans aucun doute afin de mener à bien ma future carrière de chercheur. Sa grande expertise dans le domaine du comportement des matériaux a été d'une grande aide et m'a permis de prendre du recul sur ces trois dernières années afin d'adopter une vision différente des choses. I could not forget to thank Professor Thomas D. SEWELL, from Mizzou University, who also agreed to be reporter for this thesis. It was an honour to have him as a member of the jury and to benefit from his undeniable experience in the field of energetic materials. In particular, I would like to thank him for his report, which I enjoyed reading a lot. Receiving such comments from Professor SEWELL was not something I would have imagined at the beginning of this thesis.

iii Je souhaite enfin adresser mes remerciements à Mme. Maurine MONTAGNAT, Directrice de recherche CNRS à l'Institut des Géosciences de l'Environnement pour son examen de mon travail. Ses multiples remarques m'ont fait comprendre que l'emploi de nouvelles approches dans le but de mieux comprendre le comportement de la matière représente une problématique commune à bien des domaines. Je tiens à remercier dans un premier temps M. Christophe DENOUAL, qui m'a accueilli dans son laboratoire il y a de c ¸a 3 ans et demi pour mon stage de fin d'études lors duquel j'ai effectué mes premiers pas dans ce sujet qui me paraissaît (et me paraît encore aujourd'hui) bien complexe. Sa grande culture scientifique, aussi bien en mécanique qu'en physique, son immense capacité d'écoute et surtout la confiance qu'il m'a donnée ont été essentielles et ont mené à de longues discussions toutes plus enrichissantes les unes que les autres. Il a été le premier à me présenter la notion d'approche multiéchelle qui a constitué le fil directeur de cette thèse et permis l'instauration d'un dialogue fort entre physique et mécanique des milieux continus. Enfin, je tiens à dire que je lui suis très reconnaissant de m'accueillir au sein de son laboratoire et j'espère sincèrement que iv les années à venir seront aussi fructueuses et enrichissantes que celles passées, mais je n'en doute pas.

Mes remerciements vont également à M. Jean-Bernard MAILLET qui a été très présent dès le début de mon stage et tout au long de la thèse.

Ayant une formation à dominante mécanique des milieux continus, j'ai découvert grâce à lui et ses grandes connaissances scientifiques le monde de la physique et notamment la méthode de la dynamique moléculaire qui m'était inconnue auparavant. Il a été une réelle source de motivation et ses multiples suggestions de nouvelles pistes à explorer ont selon moi été essentielles pour avancer dans ce long projet. Je lui suis aussi très reconnaissant pour m'avoir encouragé à communiquer autour de mon travail, que ce soit lors de conférences ou de publications.

Je remercie vivement M. Olivier CASTELNAU et Mme. Katell DER-RIEN, respectivement directeur et co-directrice de cette thèse pour leur suivi régulier lors de ces trois longues années. Bien que leur domaine de prédilection se situe sur des échelles de temps et d'espace bien supérieures à celles abordés dans ce travail, il a été néanmoins très intéressant de discuter avec eux et d'avoir un point de vue extérieur, permettant de prendre parfois du recul sur mes résultats et de ne pas toujours avancer "tête dans le guidon". Je les remercie également pour leurs connaissances scientifiques concernant ces échelles supérieures auxquelles j'ai tardivement commencé à m'intéresser et sur lesquelles je serai certainement amené à travailler dans un futur proche.

J'aimerais remercier tout particulièrement M. Hervé TRUMEL, Ingénieur de recherche au CEA le Ripault, que j'ai eu l'occasion de voir régulièrement au cours de ces trois années. J'ai fortement apprécié nos discussions scientifiques autour de ma thématique et je suis sincèrement convaincu que ses conseils ont fortement contribué au bon déroulement v de cette thèse.

Ces trois années passées au CEA auraient été bien différentes sans toutes les personnes avec qui j'ai pu discuter ou même travailler. J'aime- At CEA/DAM, the study of materials under large strain rate loadings is a topic of much interest. Especially, building constitutive laws for plastic bonded explosives (PBX) at the macroscopic scale is one of the main goals. It proves to be very challenging. PBX consists in the association of polycrystalline aggregates of Energetic Materials with a polymer binder.

Experiments at the macroscopic scale (polycrystalline aggregate + polymer) are well controlled, as well as microscopic characterizations with X-ray diffraction experiments (at the scale of the single crystal). However, at intermediate scale, the characterization of the single and polycrystals behavior represent a difficulty for some materials, due to security constraints as well as complex synthesis pathways (especially for TATB, see below). Energetic materials, usually found as molecular crystals, exhibit a complex mechanical behavior which has to be understood in order to build representative constitutive laws at the mesoscopic scale. Since such constitutive laws cannot be precisely obtained from macroscopic experiments, one has to understand the intermediate steps of the deformation, i.e. from the microscopic scale to the mesoscopic scale. However, experimental techniques are not often available to obtain a mechanical insight on such materials at these intermediate scales. Thus, the only way is the use of simulation techniques.

The modeling of the behavior of materials has a long and rich history.

However, during the last decades, technological improvement in electronic devices has led to huge progress in high performance computing with the appearance of many supercomputing facilities at, e.g., CEA/DAM with the TERA-1000 supercomputer. Such installations, which give the 1 possibility to perform massively parallel simulations, offer a new vision of materials modeling. Indeed, due to the increase of computing capacities, constitutive laws at the mesoscopic scale now aim at explicitly taking into account the evolution of microstructure and thus, the underlying physical phenomena that occur at the atomistic scale (see the "Exascale Co-design Center for Materials in Extreme Environments (ExMatEx)" project, under the leadership of Timothy Germann (LANL) and Jim Belak (LLNL)).

The only way to inform these models is to adopt a multiscale approach.

However, to address such a challenge, a dialogue between scales has to be set up. In this thesis, molecular dynamics (MD) is chosen to characterize the mechanichal behavior of the energetic material 1,3,5-triamino-2,4,6trinitrobenzene (TATB) at the microscopic scale. The implementation of several mechanical diagnoses in our MD code "STAMP" greatly improves our understanding of the underlying deformation mechanisms and lays the foundations for a future multiscale model, opening the way to a better understanding of very anisotropic crystalline materials. In the following, we introduce some general aspects on energetic molecular crystals before focusing on TATB. Finally, a brief outline of the work presented in this thesis is presented.

Energetic molecular crystals

In material science, molecular crystals occupy a special place. Similarly to metals or ceramics for example, they can be characterized through their crystal structure defined by the positions and orientations of molecules within the primitive lattice. Indeed, these materials usually consist in a regular arrangement of molecules, defining a periodic structure with low symmetry and complex molecular stacking. Molecular crystals find applications in multiple areas among which pharmaceutical [START_REF] Wildfong | Towards an understanding of the structurally based potential for mechanically activated disordering of small molecule organic crystals[END_REF][START_REF] Bates | Assessment of defects and amorphous structure produced in raffinose pentahydrate upon dehydration[END_REF][START_REF] Feng | Toward understanding the evolution of griseofulvin crystal structure to a mesophase after cryogenic milling[END_REF] and explosive compounds, which have been widely investigated during the past decades. For example, Lei et al. [START_REF] Lei | Mesoscale modeling of dislocations in molecular crystals[END_REF] studied the dislocations evolution in sucrose and paracetamol molecular crystals, while energetic materi-als (EM) such as cyclotrimethylene trinitriamine (RDX), pentaerythritol tetranitrate (PETN), cyclotetramethylene tetranitramine (HMX) and triamino trinitrobenzene (TATB) were subjected to a large amount of studies through Molecular Dynamics (MD) simulations [START_REF] Mathew | Molecular conformational stability in cyclotrimethylene trinitramine crystals[END_REF][START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF][START_REF] Mathew | Nanoindentation of the triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene: A molecular dynamics study[END_REF][START_REF] Kroonblawd | Theoretical determination of anisotropic thermal conductivity for crystalline 1,3,5-triamino-2,4,6trinitrobenzene (tatb)[END_REF][START_REF] Kroonblawd | Theoretical determination of anisotropic thermal conductivity for initially defect-free and defective tatb single crystals[END_REF][START_REF] Kroonblawd | Characteristics of energy exchange between inter-and intramolecular degrees of freedom in crystalline 1,3,5-triamino-2,4,6-trinitrobenzene (tatb) with implications for coarse-grained simulations of shock waves in polyatomic molecular crystals[END_REF] and nanoindentation experiments [START_REF] Gallagher | Dislocation slip systems in pentaerythritol tetranitrate (petn) and cyclotrimethylene trinitramine (rdx)[END_REF][START_REF] Elban | Vickers hardness testing of sucrose single crystals[END_REF][START_REF] Ramos | Direct observation of plasticity and quantitative hardness measurements in single crystal cyclotrimethylene trinitramine by nanoindentation[END_REF][START_REF] Taw | The mechanical properties of as-grown noncubic organic molecular crystals assessed by nanoindentation[END_REF]. These studies are mainly intended to investigate the mechanical response of complex molecular crystals, poorly understood in some cases, but of critical importance.

Crystalline molecular explosives represent the main ingredient of engineered energetic formulations. Since such materials are related to numerous applications, a very precise understanding of their behavior under a wide range of thermomechanical conditions is required. A comprehensive knowledge of the first steps of deformation at the microscopic level is thus necessary if one wants to predict the elasticity, plasticity, failure and phase transitions that link mechanical response to detonation performance. In 2015, for the 40 th anniversary issue of Propellants, Explosives, Pyrotechnics (PEP), Hooks et al. published a review of elasticity in molecular explosives [START_REF] Hooks | Elasticity of crystalline molecular explosives[END_REF]. They described the elastic response of PETN, RDX and HMX as well as comparisons between computational and experimental results and discussed the case of TATB, for which only simulation results are available. This review highlights the challenge of studying the elastic response of energetic materials for both experiments and simulations.

What happens between reversible elastic deformation (the first step of a mechanical stimuli) and macroscopic behavior (up to the transition to failure) is in part governed by irreversible plastic flow (we will demonstrate in this thesis that other mechanisms akin to phase transformation must be mentioned). The latter is triggered once the stress reaches a level called yield stress, which has become the center of all the attention during the last decades.

In an editorial note from Hooks et al. in PEP [START_REF] Hooks | Plasticity in crystalline molecular explosives -a key to unraveling unpredictable responses[END_REF], the authors recall the first answer of Marc Cawkwell (staff member at Los Alamos National Laboratory) about this topic: "What we know, really know, about plasticity in molecular explosives, could be written on a beer mat with space left over", reflecting that there is an undeniable lack of information concerning molecular explosives plasticity.

The ability of EM to remain stable during shock loading defines their sensitivity, a critical property for their safety. This sensitivity, a direct consequence of the underlying chemistry, is highly dependent on thermally activated processes. Besides, complex mechanisms such as localizations of deformation [START_REF] Davis | High explosives: The interaction of chemistry and mechanics[END_REF] may lead to the generation of hot-spots [START_REF] Mathew | Nanoindentation of the triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene: A molecular dynamics study[END_REF][START_REF] Field | Hot spot ignition mechanisms for explosives[END_REF][START_REF] Armstrong | Investigation of hot spot characteristics in energetic crystals[END_REF][START_REF] Nabarro | Highly shocked polymer bonded explosives at a nonplanar interface: Hot-spot formation leading to detonation[END_REF].

Such hot-spots are known to be triggered by local plastic deformation, dislocation pile-up [START_REF] Mathew | Nanoindentation of the triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene: A molecular dynamics study[END_REF], void collapse [START_REF] Eason | Molecular dynamics simulations of the collapse of a cylindrical pore in the energetic material α-rdx[END_REF][START_REF] Austin | Direct numerical simulation of shear localization and decomposition reactions in shock-loaded hmx crystal[END_REF] and defects nucleation [START_REF] Field | Hot spot ignition mechanisms for explosives[END_REF].

All these phenomena generally lead to a localized heating, which directly affects the material sensitivity. Irreversible deformation mechanisms thus need to be accurately described for a predictive modeling of ignition. Dislocations in metals are known to nucleate, interact with each other as well as with other defects such as voids, interfaces and grain boudaries, making crystal plasticity very complex. With molecular crystals comes further complexity concerning dislocation activity since this class of materials generally exhibits an exceptional brittleness. In addition, slip systems in molecular crystals are very difficult to identify due to their low symmetry and the possible overlap between molecules, occuring when neighboring molecular planes glide relatively to each other.

Although several experimental methods have been shown to be effective during the past [START_REF] Connick | Dislocation etching of cyclotrimethylene trinitramine crystals[END_REF][START_REF] Halfpenny | Dislocations in energetic materials: I. the crystal growth and perfection of pentaerythritol tetranitrate (petn)[END_REF][START_REF] Gallagher | Dislocation slip systems in pentaerythritol tetranitrate (petn) and cyclotrimethylene trinitramine (rdx)[END_REF][START_REF] Ramos | Direct observation of plasticity and quantitative hardness measurements in single crystal cyclotrimethylene trinitramine by nanoindentation[END_REF][START_REF] Eddleston | Transmission electron microscopy of pharmaceutical materials[END_REF], their revision and application to energetic materials should be undertaken while waiting for the emergence of new methods. Furthermore, the construction of interatomic potentials for MD simulations adapted to the study of EM has helped a lot to investigate mechanisms of the irreversible deformation in molecular crystals such as RDX [START_REF] Cawkwell | Homogeneous dislocation nucleation in cyclotrimethylene trinitramine under shock loading[END_REF][START_REF] Mathew | Slip asymmetry in the molecular crystal cyclotrimethylenetrinitramine[END_REF][START_REF] Pal | Rotational defects in cyclotrimethylene trinitramine (rdx) crystals[END_REF] and TATB [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF]. However, such studies need to be pursued in order to obtain the main elements necessary for the construction of a constitutive law at the mesoscopic scale.

Obviously, the behavior of the PBX has to be distinguished from the behavior of energetic material single and polycrystals at the grain scale. However, due to the weak properties of polymeric binders, the crystal behavior is supposed to be predominant, justifying the multiscale modeling from the microscopic to the mesoscopic scale.

An intriguing case: TATB

This thesis focuses on the investigation of the mechanical response under deformation of the energetic molecular crystal 1,3,5-triamino-2,4,6trinitrobenzene (TATB). In spite of numerous studies devoted to a better understanding of TATB behavior under moderate loadings, the mechanisms of its irreversible deformation still remain poorly understood. Due to the difficulties to obtain single crystals of size compatible with standard mechanical testing machines as well as security constraints, no experimental data is available on, e.g., single TATB crystal under pure uniaxial tension. For these reasons, the numerical modeling of TATB single crystal plays an essential part in studying its macroscopic behavior. TATB crystallizes in a triclinic cell of space group P1 [START_REF] Cady | The crystal structure of 1,3,5triamino-2,4,6-trinitrobenzene[END_REF][START_REF] Plisson | Equation of state of 1,3,5-triamino-2,4,6trinitrobenzene up to 66 gpa[END_REF] and exhibits very anisotropic thermo-mechanical and chemical behavior [START_REF] Kroonblawd | Theoretical determination of anisotropic thermal conductivity for crystalline 1,3,5-triamino-2,4,6trinitrobenzene (tatb)[END_REF][START_REF] Kroonblawd | Theoretical determination of anisotropic thermal conductivity for initially defect-free and defective tatb single crystals[END_REF] due to the arrangement in layers of the two C 6 H 6 N 6 O 6 molecules present in its unit cell (see Figure 1.1). Such anisotropic properties are also the consequences of the directional interactions within the molecular crystal, which give to the strength a high sensitivity toward the direction of loading (up to a factor 20 for shear), as indicated by computations of the fourth order stiffness tensor C [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF][START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF][START_REF] Lafourcade | Dislocation core structure at finite temperature inferred by molecular dynamics simulations for 1,3,5-triamino-2,4,6-trinitrobenzene single crystal[END_REF].

Experiments have pointed out that TATB grown single crystals are noticeably defective [START_REF] Phillips | Some observations on the structure of tatb[END_REF][START_REF] Hoffman | Density distributions in tatb prepared by various methods[END_REF][START_REF] Mallik | Understanding the solubility of triamino-trinitrobenzene in hydrous tetramethylammonium fluoride: a first principles molecular dynamics simulation study[END_REF][START_REF] Guerain | Structural defect evolution of tatb-based compounds induced by processing operations and thermal treatments[END_REF], with presence of porosities (see Fig-

ure 1.
2), twinned structures, etc. Until recently, TATB single crystals of size compatible with standardized mechanical experiments (e.g. tensile loading) could not be obtained, mainly due to the complexity of the synthesis process. However, Taw et al. [START_REF] Taw | The mechanical properties of as-grown noncubic organic molecular crystals assessed by nanoindentation[END_REF] managed to obtain quasihexagonal plate-like crystals of up to 1 mm width with a minimal presence of impurities. They were able to perform nanoindentation experiments on as-grown TATB single crystals with sizes varying between 200 µm and 600 µm which represents a major step in the experimental area for this material.

MD simulations of TATB single crystal are widely used and play a crucial part in predicting its macroscopic behavior. Kroonblawd et al. [START_REF] Kroonblawd | Theoretical determination of anisotropic thermal conductivity for crystalline 1,3,5-triamino-2,4,6trinitrobenzene (tatb)[END_REF][START_REF] Kroonblawd | Theoretical determination of anisotropic thermal conductivity for initially defect-free and defective tatb single crystals[END_REF][START_REF] Kroonblawd | Characteristics of energy exchange between inter-and intramolecular degrees of freedom in crystalline 1,3,5-triamino-2,4,6-trinitrobenzene (tatb) with implications for coarse-grained simulations of shock waves in polyatomic molecular crystals[END_REF] focused their work on the thermal behavior of TATB and relaxation of idealized hot spots. Dislocations are known to be the main mechanism that controls plasticity in crystalline materials but their existence in TATB has never been assessed. The development of predictive MD potentials [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF][START_REF] Rai | Monte carlo simulations of 1,3,5-triamino-2,4,6-trinitrobenzene (tatb): Pressure and temperature effects for the solid phase and vapor-liquid phase equilibria[END_REF] is the first step toward the prediction of dislocation core structures and their mobility. Mathew et al. computed the first generalized stacking fault surfaces of TATB [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] and presented MD simulations of its elasticplastic response under displacement-controlled nanoindentation [START_REF] Mathew | Nanoindentation of the triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene: A molecular dynamics study[END_REF].

The latter results are shown to be in very good agreement with the recent experimental work of Taw et al. [START_REF] Taw | The mechanical properties of as-grown noncubic organic molecular crystals assessed by nanoindentation[END_REF]. Still, the understanding of plasticity mechanisms in TATB single crystal is unsatisfactory compared to other explosive materials such as RDX or HMX.

Manuscript outline

The principal aim of this thesis is to build a constitutive law for the mechanical behavior of the energetic material 1,3,5-triamino-2,4,6trinitrobenzene (TATB). Because of the recent improvements in high performance computing, constitutive laws at the mesoscopic scale are now able to take into account the materials microstructure and thus the underlying physics at the atomic scale. In this work, classical MD simulations are used to provide a better understanding of the mechanical behavior of TATB single crystal. A characterization and a selection of the essential ingredients at the microscopic scale is proposed in order to supply a mesoscopic model. Furthermore, the implementation of several mechanical tools allows the dialogue between MD simulations and continuum mechanics when one to one comparison on single crystal simulations are performed. The multiscale modeling is at the core of this thesis and the focus is put on the importance of this "bottom-up" approach. An understanding of the elementary mechanisms of the irreversible deformation of materials at the atomic scale is crucial to build a constitutive law at the mesoscopic scale in which relevant ingredients from the lower scale are taken into account. Such a work requires a lot of computing ressources.

Indeed, all subsequent simulations presented in the manuscript have been performed with parallel computations, with a total cumulative time of 20 millions CPU hours.

The thesis is organized as follows. Chapter 2 is dedicated to a literature review on the studies that have been undertaken to understand the general behavior of the explosive compound 1,3,5-triamino-2,4,6trinitrobenzene. We first introduce a review of physical properties of TATB single crystal obtained by both experiments and simulations. We then introduce the reader to some physical aspects that concern the macroscopic properties of a TATB-based plastic bonded explosive (PBX).

Chapter 3 describes the simulation methods that have been employed during this thesis to bridge microscopic and mesoscopic scales which is the guiding thread of this work. We first recall the basic principles of classical MD which has proved indispensable to obtain the main elements necessary for the construction of solid foundations for a constitutive law at the mesoscopic scale. The force field used is our MD simulations is presented together with basic simulations that help for its validation such as the equation of state (EOS) and elastic properties of the single crystal. We then turn to the context of large deformations in molecular dynamics. We present a method to prescribe deformation paths to a MD simulation cell which has been implemented for this study, as well as a local deformation gradient tensor computation and its polar decomposition. These two diagnoses are typical to continuum mechanics.

Their use thus improves the dialogue between atomic and mesoscopic scales even if different simulation techniques have to be used. Finally, we present the Peierls-Nabarro-Galerkin (PNG) model (a semi-continuous method) that allows us to describe the dislocation core structure in details.

In chapter 4, we focus on the dislocation core structure at finite temperature in basal planes (noted P1 and P2) of TATB single crystal.

Investigation of gliding properties in P1 and P2 basal planes is proposed by coupling two simulation techniques. A new method is proposed to compute the generalized stacking fault surfaces (denoted γ-surfaces) as well as the elastic tensor at finite temperature through MD simulations;

these two physical properties are then used in the PNG model to study the dislocation core structures of [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] and [010] edge and screw dislocations at ambient temperature and for pressures up to 10 GPa.

In order to study the mechanisms of deformation, we present in chapter 5 a procedure, using the deformation paths prescription presented in chapter 3. Thanks to this simulation setup, a yield stress surface for TATB single crystal under pure shear is presented, leading to the identification of three distinct behaviors, each one associated with a direction of loading. The first one is related to the results of chapter 4 when pure shear is applied to the basal plane. Second, we predict the nucleation of dislocations in non-basal planes associated with an unusual local dilatancy process. And third, every deformation containing a basal-plane compressive component involves a buckling-like deformation. We then relate it to a deformation path that allows for a perfect twinning of the TATB triclinic cell, found through the computation of a twinning energy landscape analogous to a γ-surface.

The gathering of all the physical properties calculated via MD simulations led us to a better understanding of the irreversible deformation phenomena in TATB single crystal. Chapter 6 is dedicated to a discussion and results on the modeling of TATB single and polycrystal at the mesoscopic scale using an element free Galerkin method (EFG). We show that for materials with large anisotropy, such as TATB, two equations of state (EOS) should be used, depending on the loading type (by prescribed pressure, or deformation).

Finally, chapter 7 is intended to bring together the key elements and conclusions of this work and to propose some perspectives for future work. [START_REF] Mathew | Nanoindentation of the triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene: A molecular dynamics study[END_REF][START_REF] Field | Hot spot ignition mechanisms for explosives[END_REF][START_REF] Armstrong | Investigation of hot spot characteristics in energetic crystals[END_REF][START_REF] Nabarro | Highly shocked polymer bonded explosives at a nonplanar interface: Hot-spot formation leading to detonation[END_REF]. Ces points chauds sont connus pour être déclenchés par la déformation plastique locale, l'empilement de dislocations [START_REF] Mathew | Nanoindentation of the triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene: A molecular dynamics study[END_REF], la fermeture de pores [START_REF] Eason | Molecular dynamics simulations of the collapse of a cylindrical pore in the energetic material α-rdx[END_REF][START_REF] Austin | Direct numerical simulation of shear localization and decomposition reactions in shock-loaded hmx crystal[END_REF] ou encore à la nucléation de défauts [START_REF] Field | Hot spot ignition mechanisms for explosives[END_REF]. [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] et ont présenté des simulations de DM de la réponse élasto-plastique sous nanoindentation du TATB [START_REF] Mathew | Nanoindentation of the triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene: A molecular dynamics study[END_REF].

Introduction en franc ¸ais

Ces derniers résultats s'avèrent être en très bon accord avec le récent travail expérimental de Taw et al. [START_REF] Taw | The mechanical properties of as-grown noncubic organic molecular crystals assessed by nanoindentation[END_REF].

Plan du manuscrit

Le but de cette thèse est de mettre en place une loi de comportement pour la réponse mécanique du matériau énergétique 1,3,5-triamino-

CHAPTER 2 REVIEW OF THE SCIENTIFIC LITERATURE

This chapter is dedicated to a short review of existent work on the TATB behavior. We first recall various studies that have been conducted on the single crystal in section 2.1 with a focus on the few available experimental results and a presentation of atomistic models developed for TATB with related results on thermal and mechanical properties obtained through MD simulations. We then briefly present in section 2.2 both experimental and computational results on the macroscopic behavior of TATB-based plastic-bonded explosives (PBX). on a 200 µm sample. They were the first to determine that TATB crystallizes in a triclinic cell of space group P1, with a graphitic-like molecular stacking. They identified the strong in-plane hydrogen-bond network that is believed to provide TATB its high structural stability as long as the weak van der Waals interactions between molecular layers, leading to strong anisotropic properties for this molecular crystal. A large amount of work has been undertaken to study the structural properties of TATB and the potential existence of other polymorphs [START_REF] Kolb | Growth of 1,3,5-triamino-2,4,6trinitrobenzene (tatb) i. anisotropic thermal expansion[END_REF][START_REF] Son | Dynamic observation of a thermally activated structure change in 1,3,5-triamino-2,4,6-trinitrobenzene (tatb) by second harmonic generation[END_REF][START_REF] Filippini | The crystal structure of 1,3,5-triammo-2,4,6-trimtrobenzene: Centrosymmetric or noncentrosymmetric?[END_REF][START_REF] Davidson | Stubborn triaminotrinitrobenzene: Unusually high chemical stability of a molecular solid to 150 gpa[END_REF][START_REF] Phillips | Some observations on the structure of tatb[END_REF] in addition to the one found by Cady and Larson. However, no other forms than the triclinic structure have been identified ever since. In 2007, Stevens et al. [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] determined a hydrostatic compression curve for TATB up to

Structure and properties of TATB single crystal

13 GPa through powder X-ray diffraction in conjunction with a diamond anvil cell. Similarly, in 2017, Plisson et al. [START_REF] Plisson | Equation of state of 1,3,5-triamino-2,4,6trinitrobenzene up to 66 gpa[END_REF] obtained the equation of state of TATB for pressures between 0 and 66 GPa using synchrotron X-ray diffraction experiments in a diamond anvil cell. The experimental unit cells parameters are reported in Table 2.1. Additionally, spectroscopic studies (both experimental and theoretical) have shown the strenghtening of the in-plane hydrogen-bond network with pressure [START_REF] Pravica | Infrared study of 1,3,5-triamino-2,4,6-trinitrobenzene under high pressure[END_REF][START_REF] Pravica | High-pressure far-and mid-infrared study of 1,3,5-triamino-2,4,6trinitrobenzene[END_REF][START_REF] Manaa | Nearly equivalent inter-and intramolecular hydrogen bonding in 1,3,5-triamino-2,4,6trinitrobenzene at high pressure[END_REF].

Atomistic potentials for TATB

Various force fields (FFs) have been proposed for the modeling of TATB behavior through MD simulations as well as modifications and post-parametrizations [START_REF] Gee | Ab initio based force field and molecular dynamics simulations of crystalline tatb[END_REF][START_REF] Rai | Monte carlo simulations of 1,3,5-triamino-2,4,6-trinitrobenzene (tatb): Pressure and temperature effects for the solid phase and vapor-liquid phase equilibria[END_REF][START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF][START_REF] Kroonblawd | Theoretical determination of anisotropic thermal conductivity for crystalline 1,3,5-triamino-2,4,6trinitrobenzene (tatb)[END_REF]. These FFs include reactive, non-reactive properties and some of them treat TATB molecules as rigid. Therefore, when considering rigid molecules, the ability to treat conformational flexibility is removed as for example with the Rai [START_REF] Rai | Monte carlo simulations of 1,3,5-triamino-2,4,6-trinitrobenzene (tatb): Pressure and temperature effects for the solid phase and vapor-liquid phase equilibria[END_REF] FF which has been used in the past decade to compute properties such at thermal expansion coefficients and elastic constants. A brief review of these studies is proposed in the following.

Mechanical properties

The anisotropic graphitic-like molecular stacking of TATB single crystals is suspected to be responsible of its very anisotropic mechanical properties. However, due to the difficulties to obtain perfect single crystals of size compatible with standard mechanical testing techniques as well as security constraints, the majority of mechanical properties for TATB single crystal have been predicted through numerical simulations, using Density Functional Theory (DFT), Monte-Carlo simulations (MC) and classical Molecular Dynamics (MD).

No experimental determination of TATB elastic constants have been realized so far. Bedrov et al. [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF] computed the isothermal fourth order elastic tensor for TATB by using the Parrinello-Rahman strain-strain fluctuation formula [START_REF] Parrinello | Strain fluctuations and elastic constants[END_REF] (described in appendix A.2) for pressures up to 10 GPa on the 300 K isotherm, following a Molecular Dynamics-Monte Carlo (MD-MC) approach [START_REF] Borodin | Development of quantum chemistry-based force fields for poly(ethylene oxide) with many-body polarization interactions[END_REF][START_REF] Sewell | A molecular dynamics simulation study of elastic properties of hmx[END_REF]. The elastic anisotropy of TATB is preserved with pressure with a large increase of C 11 and C 22 compared to C 33 , with a similar evolution for C 66 in comparison to C 44 and C 55 . They also computed the temperature dependency of the stiffness tensor for temperatures between 200 K and 400 K at zero pressure. The effect of temperature on elastic constants is a global softening on compressive and dilational elements while shear components are found to be nearly independent on temperature.

Another set of elastic constants for TATB single crystal was obtained by Mathew et al. [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] at ambient (T,P) conditions using the same approach but with a modified version [START_REF] Kroonblawd | Theoretical determination of anisotropic thermal conductivity for crystalline 1,3,5-triamino-2,4,6trinitrobenzene (tatb)[END_REF] of the potential used by nisms of the deformation in TATB single crystal is still limited. The layered structure gives TATB some similarities with graphite, boron nitride and other layered materials. When it comes to plasticity in materials of any kind, dislocations motion is the first mechanism to be considered. For layered materials, as discussed by Delavignette [START_REF] Delavignette | Dissociation and plasticity of layer crystals[END_REF], plastic deformation is essentially governed by dislocations glide along planes parallel to the layers (also called basal planes). In the majority of layered crystals, binding energies between atoms/molecules within the layers (strong hydrogen bonding for TATB) are way stronger than interactions between successive layers (van der Waals type for TATB). While some exceptions exist, dislocations in layered crystals are subjected to a dissociation phenomena in two or more partials [START_REF] Delavignette | Dissociation and plasticity of layer crystals[END_REF][START_REF] Amelinckx | Dislocations in particular structures[END_REF]. This dissociation is known to be controlled by the stacking fault energies (γ SF ) along the gliding direction. In order to study these properties for TATB single crystal, Mathew et al. [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] computed the generalized stacking fault energies surfaces (γ-surfaces) for the (001) basal plane. Due to the molecular arrangement in TATB, the glide in two basal planes was investigated. They reported low-energy stacking faults at 0 GPa indicating easy glide of dislocations in the basal plane as well as low critical stresses for dislocation motion. Furthermore, Zhang investigated the sliding properties of a TATB single layer through density functional theory [START_REF] Zhang | Investigation of the slide of the single layer of the 1,3,5-triamono-2,4,6-trinitrobenzene crystal: Sliding potential and orientation[END_REF].

Thermal properties

Several thermal properties of the single crystal were investigated in the past decade. Kroonblawd and Sewell extended an existing rigid-bond force field [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF] to include bond stretching and three-center angle bending potentials for simulations requiring fully flexible molecules [START_REF] Kroonblawd | Theoretical determination of anisotropic thermal conductivity for crystalline 1,3,5-triamino-2,4,6trinitrobenzene (tatb)[END_REF]. They predicted the thermal conductivity of crystalline TATB along three directions, resulting in a significant anisotropy with a thermal conductivity within the molecular layers approximately 1.7 times greater than along the interplanar direction, which differs significantly from values observed in others crystalline explosives such as RDX and HMX. Later, this work was extended to investigate finite-size and molecular vacancy effects on the anisotropic thermal conductivity [START_REF] Kroonblawd | Theoretical determination of anisotropic thermal conductivity for initially defect-free and defective tatb single crystals[END_REF]. Conduction within the layers is invariant with respect to the size of simulation cell in these directions while in the transverse direction the conduction is found to decrease with increasing size of supercell. Introduction of randomly distributed vacancies in crystalline TATB led to a decrease of the conduction for directions normal to the molecular layers. Finally the authors studied the sensitivity of predictions for anisotropic thermal conductivity λ to intramolecular terms of the force field. They also investigated the P = 0.0 GPa isobaric temperature dependence of λ over the temperature interval 200 K ≤ T ≤ 700 K as well as the T = 300 K isothermal pressure dependence of λ over the pressure interval 0.0 GPa ≤ P ≤ 2.5 GPa [Kroonblawd, M.P. and Sewell, T.D.]. A decrease in 1/T of λ is predicted with increasing temperature with a maximum anisotropic factor of 0.7 at 200 K while λ is predicted to increase linearly with increasing pressure.

More recently, Kroonblawd and Sewell investigated the anisotropic hot spot relaxation in crystalline TATB on length scales up to ≈ 10 nm using MD simulations with comparison with analytic solutions to the continuum-level diffusive heat equation [START_REF] Kroonblawd | Anisotropic relaxation of idealized hot spots in crystalline 1,3,5-triamino-2,4,6-trinitrobenzene (tatb)[END_REF]. They considered two con-ductions directions, according to the anisotropic thermal conductivity observed in their previous work, i.e. within and between the molecular layers. Excellent agreement between MD calculations and analytical solution were found for the conduction normal to the basal plane while some discrepancies for the other directions were smoothed by introducing a temperature dependence of thermal conductivity.

Other studies related to thermal and liquid-state properties of TATB were conducted. Kroonblawd et al. investigated the rate of energy exchange between internal and external degrees of freedom for the crystalline explosive using MD simulations with flexible molecules, leading to the adjustment of parameters for coarse-grained dissipative particles dynamics (DPDE) model for TATB [START_REF] Kroonblawd | Characteristics of energy exchange between inter-and intramolecular degrees of freedom in crystalline 1,3,5-triamino-2,4,6-trinitrobenzene (tatb) with implications for coarse-grained simulations of shock waves in polyatomic molecular crystals[END_REF]. The melt curve and liquid-state transport properties of TATB were predicted as functions of pressure and temperature by Mathew et al. [START_REF] Mathew | Predicted melt curve and liquid-state transport properties of tatb from molecular dynamics simulations[END_REF] using MD simulations with a nonreactive force field with flexible molecules.

Macroscopic properties of a TATB-based PBX

Microstructure

At the macroscopic scale, TATB grains (polycrystalline aggregates) are combined with a polymeric binder, whose presence is less than 5% in volume. The manufacturing process is done in two steps: first, TATB yellow powder is coated with the liquid-phase binder while being granulated in millimeter-beads. After a drying step, an isostatic pressing under vacuum in an oil bath is applied to the granular material at moderated temperature and high pressure. This results in a low-porosity PBX (less than 2.5%), with good mechanical strength, easy to process, with an isotropic microstructure [START_REF] Maiti | Irreversible volume growth in polymer-bonded powder systems: Effects of crystalline anisotropy, particle size distribution, and binder strength[END_REF][START_REF] Ambos | Numerical modeling of the thermal expansion of an energetic material[END_REF]. The final microstructure, with grains of typical size of 50 µm is represented in Figure 2.1. Residual porosity can be identified as small black points clusters, while other "mechanical" defects are very easy to detect: intragranular fractures are present within some grains (red ellipses) while straight bands (green circle), missing before the isostatic pressing, have been triggered and are thought to be the results of plastic deformation.

Mechanical behavior

The quasi-static behavior of this TATB-based material has been deeply characterized with a large amount and varied experiments: simple, monotonous and cyclic compression tests, low and high pressure tri-axial compression tests and cyclic hydrostatic compression tests have been performed [START_REF] Ambos | Numerical modeling of the thermal expansion of an energetic material[END_REF].

Under simple traction, its brittle behavior [START_REF] Blumenthal | Influence of temperature and strain rate on the mechanical behavior of pbx 9502 and kel-f 800[END_REF][START_REF] Thompson | Timetemperature analysis, tension and compression in pbxs[END_REF] leads to a rupture within the elastic domain whereas for simple compression, damage progresses slowly until rupture is provoked at 2% of longitudinal deformation, leading to a usual factor 5 on yield stresses between uniaxial traction and compression.

An evolution with pressure of the strain-level at failure is also observed. A longitudinal deformation of 1.5% for a simple compression leads to a 30 MPa yield stress while under a confining pressure of 200 MPa, the yield stress approaches 85 MPa for a longitudinal strain around 30%. This is typical of a brittle-ductile transition found in quasi-brittle materials under pressure.

Thermoelastic behavior

Young's and transverse moduli have been estimated from experiments [START_REF] Ambos | Numerical modeling of the thermal expansion of an energetic material[END_REF] with values of 7.1 GPa and 21.2 GPa, respectively, leading to a

Poisson's coefficient around 0.335. This allows the estimation of the bulk modulus (7.2 GPa) and shear modulus (2.7 GPa). Furthermore, a volumetric thermal expansion coefficient of 1.62 × 10 -4 K -1 has been found

and is comparable to the half of that of TATB single crystal presented in section 2.1.4.

A lot of work is currently underway to simulate the PBXs macroscopic behavior within a wide range of temperature and pressure. Ambos et al. [START_REF] Ambos | Numerical modeling of the thermal expansion of an energetic material[END_REF] proposed a numerical modelling of the thermal expansion of a TATB-based pressed explosive while, Gasnier et al. [START_REF] Gasnier | A fourier-based numerical homogenization tool for an explosive material[END_REF] introduced a FFT homogenization tool to compute thermoelastic properties of the same material. More recently, Gasnier [START_REF] Gasnier | Etude du comportement thermo-mécanique et de l'endommagement d'un matériau énergétique granulaire par méthodes de fourier[END_REF] started to investigate the impact of microfissuration on its macroscopic properties. However, it has been pointed out that these models do not include any of the irreversible mechanisms of the deformation of TATB single crystal since no data (neither experimental nor numerical) are available on the subject.

Conclusion

The literature review presented above has shown the difficulty and the challenge to obtain data on the mechanical behavior of TATB from experiments. Indeed, the TATB single crystal structure is known from X-ray diffraction but no mechanical testings at the mesoscopic scale have been performed so far except recently by Taw et al. [START_REF] Taw | The mechanical properties of as-grown noncubic organic molecular crystals assessed by nanoindentation[END_REF] with nanoindentation experiments. A few numerical investigations have been done by Mathew et al. who studied the gliding properties in TATB basal planes [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] and its response to nanoindentation [START_REF] Mathew | Nanoindentation of the triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene: A molecular dynamics study[END_REF]. The very first aim of this thesis is thus the investigation of the possible deformation mechanisms of TATB single crystal in order to build solid foundations for a future constitutive law at larger time and space scales.

Résumé substantiel en franc ¸ais

Lagangian and Hamiltonian formulations

The classical Molecular Dynamics (MD) method consists in studying the evolution of a physical system composed of many particles (atoms or molecules) in order to investigate its microscopic behavior under various temperature and pressure conditions [START_REF] Allen | Computer Simulation of Liquids[END_REF][START_REF] Frenkel | Understanding Molecular Simulation[END_REF]. A general basis for MD simulations can be obtained through the Lagrangian formulation of classical mechanics. The Hamilton's variational principle states that a mechanical system of Lagrangian L will follow a trajectory in the phase-space for which the integral Ldt is an extremum. Let consider N particles interacting with each other with independent generalized coordinates and velocities {q i , qi }, describing the evolution of a conservative system of which forces derive from a potential U . Then, L = L({q i }, { qi }, t) satisfies the Lagrange equations:

d dt ∂L ∂ qi - ∂L ∂q i = 0, i = 1, ..., N, (3.1)
which constitute a basis for Newton's second law. Considering q i as a cartesian coordinates component of an atom, associated with a mass m, the Lagrangian reads:

L = 1 2 m i q2 i -U ({q i }), (3.2)
leading to the change of equation 3.1 into

mq i = - ∂U ∂q i = F i , (3.3)
where F i represents the force derivating from the potential function U .

Generaly, the Lagrange function of a given system can be written as:

L = K(q, q) -U (q), (3.4)
with K and U representing kinetic and potential energy, respectively. One can also introduce the system's Hamiltonian H:

H = i qi • p i -L = K + V, (3.5)
where the generalized impulsion p i is defined by

p i = ∂L ∂ qi
. This formulation leads to the two following first-order equations of motion, known as the Hamilton's equations:

qi = ∂H ∂p i , ṗi = - ∂H ∂q i . (3.6)
In general, these equations can not be resolved in an analytical way and are discretized in order to be integrated on short time intervals ∆t of the order of femtosecond. When interactions between atoms are described through an interatomic potential U (q 1 , q 2 , ..., q i , ..., q N ), each particle undergoes the behavior of its neighbors at a time t and the force F i acting on the i th atom is computed as:

F i = -∇U , (3.7)
where ∇ is the gradient operator.

Working with rigid molecules

Depending on the chosen time step to perfom MD simulations (≃ 1.0 fs in our case), some intramolecular interactions (as bond vibrations effects for example) have typical times below intermolecular interactions ones and may be considered independently, allowing the rigid molecules framework to simulate the system. Furthermore, certain types of molecules present quasi-rigid bond properties with no possible conformational configurations. Therefore, trajectories are computed from the molecular centers of mass (COM) and orientations using the rigid body dynamics. The orientation of a rigid molecule can be fully interpretated through the quaternion (see appendix A.1), a notion first introduced by

Evans [START_REF] Evans | On the representatation of orientation space[END_REF]. In the present case, the Lagrange function reads:

L = 1 2 mV 2 + 1 2 Iω 2 -U , (3.8)
where V represents the COM velocity, m the mass, I the inertia tensor, ω the angular velocity around an axis passing through the COM and U the potential energy.

Numerical intergration of equations of motion

In order to compute the forces acting on each particle, potentials are usually truncated so that a particle only knows its neighbors located at a distance smaller than r cut called the cut-off radius. Some corrections are generally applied in order to ensure the forces and energy continuity since the energy rarely vanishes at this exact distance. During the trajectory, some methods have to be implemented to allow the search for neighbors.

The Verlet list [START_REF] Verlet | Computer experiments on classical fluids. ii. equilibrium correlation functions[END_REF] consists in storing all the neighbors that are located in a surrounding volume with a radius r v > r cut for each particle. Choosing r v large enough implies that all the neighbors at a distance smaller than r cut are contained in the Verlet list, even after the integration of the equations of motion during one time step. Furthermore, as long as the displacement of a particle remains smaller than r v -r cut , the Verlet list is unchanged.

Since all the forces acting between particles are now computed, Newton's equations of motion can be integrated. This is done by computing positions and orientations of molecules at time t + ∆t, ∆t being the timestep for the discretization of equations of motion. To do so, a Taylor expansion of particle position around time t is done:

r(t + ∆t) = r(t) + ∆t ṙ(t) + ∆t 2 2! r(t) + ∆t 3 3! ... r (t) + ... (3.9)
where r(t) is the vector position of a given particle at time t. With v(t) = ṙ(t) and a(t) = r(t) = f (t) m , with m the particle mass, this expression similarly reads:

r(t + ∆t) = r(t) + v(t)∆t + f (t) 2m ∆t 2 + ∆t 3 3! ... r (t) + ... (3.10)
To solve this equation, the Velocity-Verlet [START_REF] Verlet | Computer experiments on classical fluids. ii. equilibrium correlation functions[END_REF] algorithm was used, for precision and stability purposes. It consists in evaluating first the velocities at t + ∆t 2 before computing r and v at t + ∆t as follows:

1. v(t + 1 2 ∆t) = v(t) + 1 2 a(t)∆t 2. r(t + ∆t) = r(t) + v(t + 1 2 ∆t)∆t 3. v(t + ∆t) = v(t + 1 2 ∆t) + 1 2 a(t + ∆t)∆t (3.11)
This integrator ensures time reversibility (i.e. the inversion of the timestates sequence keeps the dynamics of the system well-defined) and is a symplectic integrator.

Thermodynamics ensembles

Studying atomic trajectories of materials through molecular dynamics simulations can help obtaining some of their macroscopic properties.

These properties are calculated by using a thermodynamic ensemble, according to the type of properties to be calculated. We briefly focus on the ensembles that have been used during this work, i.e. the microcanonical ensemble and the canonical ensembles, denoted NVE and NVT ensembles, respectively [START_REF] Tuckerman | Statistical Mechanics: Theory and Molecular Simulation[END_REF][START_REF] Diu | Physique Statistique[END_REF]. An other technique, the Parrinello-Rahman method, that includes a modification of the Lagrange equations of motion has also been used and is described in detail in the appendix A.2.

Microcanonical ensemble

The microcanonical ensemble consists in considering an isolated system with N particles, in a volume V , of energy E. By definition, isolated means that the system does not exchange neither particles nor energy with its surrounding environment. Therefore, the simulations performed in this ensemble are done at constant number of particles N , volume V and energy E, by solving Newton's equations without any temperature or pressure control. Distribution and probability of possibles states are described by a partition function Q NVE , defined as:

Q NVE = 1 N !h 3N δ E(q, p) -E 0 dqdp, (3.12)
where h is the Planck constant, dqdp an integral over the phases space, E(r, p) the configuration energy, E 0 the desired energy and δ the Kronecker symbol.

Canonical ensemble

Sometimes it is more convenient to use the NVT ensemble in which the temperature T is kept constant instead of the energy E. The canonical ensemble concerns a system exchanging energy with a thermostat that imposes a temperature T . The partition function for the canonical ensemble is not uniform and is known as the Boltzmann's distribution which is defined as:

Q NVT = 1 N !h 3N e -E(q,p) k b T dqdp, (3.13)
where k B is the Boltzmann constant and E is the system energy. In the N V T ensemble, the energy E varies with the latter measure and the temperature is maintained constant via the introduction of a thermostat. The most common methods used to introduce a thermostat are the algorithms of Berendsen et al. [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF], Nosé-Hoover [START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF] and Langevin [START_REF] Langevin | Sur la théorie du mouvement brownien[END_REF].

In the present work, all MD simulations were performed using the STAMP molecular dynamics software [START_REF] Soulard | Molecular dynamics study of the micro-spallation[END_REF]. Simulations were mainly per-formed on the supercomputer TERA1000 at CEA. Visual representation and analysis of simulations were performed using Visual Molecular Dynamics Software, ParaView and Blender.

Force field for TATB with rigid molecules

A force-field for fully flexible molecules initially developed by Borodin et al. [START_REF] Borodin | Polarizable and nonpolarizable force fields for alkyl nitrates[END_REF] (Atomistic Polarizable Potential for Liquids, Electrolytes, and Polymers) was adapted by Bedrov [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF] to model TATB. The total potential energy for this force field reads:

U tot = U NB (r) + bonds U bond (r ij ) + bends U bend (θ ijk ) + dihedrals U dihedral (φ ijkl ) + improper dihedrals U imp (φ imp ijkl ) (3.14)
where U NB is the potential part for non-bonded interactions and the sums concerning intramolecular terms are computed for covalent bonds, bends, dihedrals and improper dihedrals configurations in the system. The nonbonded potential energy consists in the sum of a two-bodies repulsiondispersion term, coulombian interactions and a polarization energy due to the interactions of induced dipoles. It is defined as:

U NB (r) =U RD (r) + U coul (r) + U pol (r) = i>j A αβ e (-B αβ r ij ) -C αβ r 6 ij + D 12 B αβ r ij 12 + i>j q i q j 4πǫ 0 r ij - 1 2 i µ i E 0 i . (3.15)
where

µ i = α i E tot i
is the induced dipole at force center i with α i the isotropic atomic polarizability, E tot i is the total electrostatic field at atomic site i uniquely due to the permanent charges q j and induced dipoles µ j , ǫ 0 is the vacuum dielectric permittivity and E 0 i is the electric field due to fixed charges only. A αβ and B αβ are the repulsion parameters while C αβ is the dispersion parameter between atoms i and j of type α and β. The term D(12/(B αβr ij )) 12 , with D = 5 × 10 -5 kcal mol -1 for any pair interaction, vanishes for typical atomic separations but becomes predominant for r ij < 1 Å, ensuring that U RD (r) is repulsive for short distances. For heteroatomic interactions, Waldman-Hagler [START_REF] Waldman | New combining rules for rare gas van der waals parameters[END_REF] combining rules are used with:

A ij = A ii A jj B 6 ij B 3 ii B 3 jj ; B ij = 2 B ii -6 + B jj -6 1 6 ; C ij = C ii C jj . (3.16)
Two versions for TATB force field were proposed by Bedrov [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF], one polarizable with no intramolecular nonbonded interactions and a non polarizable version where atomic polarizabilities were not considered.

Atomic polarizabilities were determined by fitting them to the molecular polarizability of gas-phase molecule. Partial atomic charges are computed to reproduce quantum chemistry (QC) results for the electrostatic potential on a regular grid around a gas-phase molecule. Bond lengths and equilibrium bending angles are ajusted to reproduce the molecule geometry in gas-phase from QC while bending force constants are taken from previously developed force fields. Finally, torsional and out-of-plane deformations parameters are determined to fit conformational energies of gas-phase molecule obtained from QC computations. Original parameters for the TATB force field are presented in Table 3.1. More detail on the force field parametrization is available in the original paper [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF].

In this work, we use the non-polarizable version of the force field developed by Bedrov. He indicated that when using the non-polarizable force field, empirical adjustments could be made, on partial atomic charges in particular. He found that the nonpolarizable force field with 20% scaled-charges showed the best agreement with experiment. We chose to perform all our MD simulations with this version along with the use of planar, symmetric rigid molecules, an important choice that is discussed in the next section. Thus, the parameters of Table 3 and torsion properties were not used. Instead of the Ewald summation for the computation of long-range electrostatic interactions, the Reaction Field [START_REF] Onsager | Electric moments of molecules in liquids[END_REF][START_REF] Barker | Monte carlo studies of the dielectric properties of water-like models[END_REF][START_REF] Watts | Monte carlo studies of liquid water[END_REF] approximation with a 13 Å cutoff radius and 100.0

K B θ 0 Bond type (constrained, Å) U B = 0.5K B (θ -θ 0 ) 2 kcal mol -2 rad -1 deg C-C 1.440 C-C(NH 2 )-C 144 118.4 C-N(O 2 ) 1.427 C-C(NO 2 )-C 144 121.6 C-N(H 2 ) 1.322 C-C-N(H 2 ) 140 121.7 N-O 1.230 C-C-N(O 2 ) 140 119.2 N-H 1.010 C-N-O 140 116.8 C-N-H 110 112.0 O-N-O 95 117.3 H-N-H 125 119.0 Torsion type K T,1 K T,2 Improper dihedral type K ID U T = 0.5K T,i [1 -cos(iφ)] (kcal mol -1 ) (kcal mol -1 ) U ID = 0.5K ID φ 2 (kcal mol -2 rad) C-C-C-C 0.0 5.35 C-C-C-N * 36.5 C-C-C-N 0.0 20.0 O-N-O-C * 89.3 N-C-C-N 0.0 20.0 H-N-H-C * 2.1 H-N-C-C 0.0 8.88 O-N-C-C 0.0 1.60
for the dielectric constant ε was used. A Langevin thermostat was used with a damping constant set to 1.0 ps. The time step for integrating the equations of motion was set to 1.0 fs for both isochoric-isoergic (NVE)

and isochoric-isothermal (NVT) simulations, small enough regarding the fact that simulations are performed with rigid molecules. Repulsion, dispersion and short-range electrostatic interactions were computed with a cut-off distance of 13 Å. All intramolecular non-bonded interactions were excluded in the implementation of the force field. Finally, the rigid body dynamics is computed by discretization of equations of motion using a Velocity-Verlet scheme [START_REF] Verlet | Computer experiments on classical fluids. ii. equilibrium correlation functions[END_REF].

Discussion on the choice of planar and rigid molecules

In a recent work [START_REF] Pal | Contribution of molecular flexibility to the elastic-plastic properties of molecular crystal α-rdx[END_REF], the contribution of molecular flexibility to the elastic-plastic properties of molecular crystal α-RDX has been studied.

The authors evaluate crystal-scale parameters such as elastic constants, lattice parameters, thermal expansion coefficients, stacking fault energies and Peierls-stress for the motion of a dislocation. They showed that the mechanical properties of α-RDX can be strongly affected by the flexibility of the constituent molecules. In our case, the benzene ring of TATB molecules provides the TATB molecules a certain rigidity and a planar geometry. Since TATB molecules do not exhibit consequent flexible behavior compared to the RDX ones, we have chosen to perform our MD simulations with perfectly planar, symmetric and rigid molecules. This hypothesis is justified when performing MD simulations in the NVE and NVT ensembles, however, since the intramolecular degrees of freedom are not taken in account, shock compression simulations would give incorrect results and flexiblity should be considered instead.

Equation of state

Molecule orientations and center of mass (COM) positions were obtained through Parrinello-Rahman [START_REF] Parrinello | Crystal structure and pair potentials: A molecular-dynamics study[END_REF][START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] simulations for different temperatures and pressures (see appendix A.2 for a detailed explanation of the method). P and T coupling constants were set to 0.1 ps and 1.0 fs, respectively. Starting with the experimental structure [START_REF] Cady | The crystal structure of 1,3,5triamino-2,4,6-trinitrobenzene[END_REF] to initiate Parrinello-Rahman simulations, a 3D-periodic 5a × 5b × 7c simulation cell containing 350 molecules was created. TATB cell parameters for each (T , P ) state were obtained by averaging the values over the final 150 ps of a 250 ps Parrinello-Rahman simulation. The average lattice parameters are given in Table 3.2. Uncertainties correspond to one standard deviation of the mean and were evaluated following the block averages method [START_REF] Frenkel | Understanding Molecular Simulation[END_REF].

It is to be noted that angle α, β and lattice parameter c deviate from experiment [START_REF] Cady | The crystal structure of 1,3,5triamino-2,4,6-trinitrobenzene[END_REF] for a few percent, probably due to the planeness and rigidity imposed to the molecules as suggested by Pal and Picu recently [START_REF] Pal | Contribution of molecular flexibility to the elastic-plastic properties of molecular crystal α-rdx[END_REF].

In Figure 3.1 is represented the cell lengths and cell angles evolution with P as well as the equation of state (EOS) of TATB. These results Experiment [START_REF] Cady | The crystal structure of 1,3,5triamino-2,4,6-trinitrobenzene[END_REF] [START_REF] Pineau | Cold curve of tatb[END_REF] for the EOS. Our results are found to be in a good agreement with the literature even though little differences can be observed on the β angle which rapidly evolves under 90 • contrarly to experiments and on the c length which is 4 % smaller, probably due to the planar symmetry imposed to the TATB molecules. Nevertheless, it can be observed that the agreement between simulations and experiments is improved for increasing pressures.

Quasi-linear coefficients of thermal expansion (CTE) with values of 1.95 × 10 -4 Å/K, 1.94 × 10 -4 Å/K and 3.60 × 10 -4 Å/K were deduced from simulation results for a, b and c directions, respectively. These values are consistent with the layered crystal structure but are shown to be quite different from those obtained by Bedrov [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF], particularly for the one in the c direction. Yet, experimental data [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] for CTE allow us to observe that our simulation underestimates the values for the b and c directions while it overestimates the one in the a direction. In their recent work, Pal

and Picu [START_REF] Pal | Contribution of molecular flexibility to the elastic-plastic properties of molecular crystal α-rdx[END_REF] demonstrate that the CTE can strongly depend on molecular structure and flexibility. In our case, the approximation on the planeness of TATB molecules and their rigidity (by constraining dihedral and angular flexibility of the ring, amino and nitro groups) may have significantly affected the CTE. The principal aim of this work is to analyze the mecha-nisms of deformation at constant T , and the rigid molecule hypothesis is believed to be a reasonable tradeoff between accuracy and numerical efficiency (10 times faster). Thus, the small differences in the crystal structure compared to the ones published in the literature [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF][START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF][START_REF] Cady | The crystal structure of 1,3,5triamino-2,4,6-trinitrobenzene[END_REF] may not have a strong influence on the general conclusions on the dislocation core structures of TATB, regarding the similarities between elastic constants and γ-surfaces obtained in this work compared to the precedent results.

Elastic properties

The complete fourth-order stiffness tensor C was obtained through a large range of T and P . The components were calculated by fitting a second-order polynomial on a set of elastic-energy vs. strain equations :

E(V, ε) = E(V 0 , 0) + V 0 3 i,j=1 (σ ij ) 0 ε ij + V 0 2 3 i,j,k,l=1 C ijkl ε ij ε kl + ... (3.17)
where Hooke's law expression is contained in the third term, expressing the second-order Cauchy stress tensor σ as a function of the linearized strains tensor ε and C, with V 0 the cell volume. These elastic-energy vs. strain curves were obtained through MD simulations by successive homogeneous deformations of a 3D-periodic 20a × 20b × 20c simulation cell. No relaxation was performed after deformation, and therefore the instantaneous (or adiabatic) elastic response is measured. This measure is however comparable to the relaxed (or isothermal) response since the prescribed deformations are too small to induce an important atomic rearrangement, with the noticeable advantage to get rid of the noise that considerably hinders the measure of an energy variation. Moreover, a system of 16 000 molecules was taken to improve the ensemble average and thus the accuracy of the measure.

In equation (3.17), E stands for potential energy. The system is con-sidered to be isentropic since no relaxation at all is performed after the strain is applied, i.e. the atomic positions remains unchanged. What is computed is consequently the isentropic elastic constants C S ij , i.e., the stiffness implied in acoustic waves for large wavelengths (see Wallace [START_REF] Wallace | Statistical Physics Of Crystals And Liquids[END_REF]).

For practical reasons, the Voigt notation C ij is used in the remaining of this paper. Two errors have to be considered here. The first comes from the average of multiple elastic-energy vs. strain curves (black dots) and is lower than 2 % while the second is due to the fitting (continuous lines) and is always lower than 1.2 %.

Results for 300 K and 0 GPa

The complete fourth-order elastic tensor at ambient conditions, which clearly reflects the anisotropic behavior of TATB, is given in Table 3.3.

This anisotropy is consistent with the crystal structure since deformations along a and b directions involve strong hydrogen bonding while only weak formula [START_REF] Parrinello | Strain fluctuations and elastic constants[END_REF] (see Table 2.2) which consists in analyzing the cell parameters fluctuations over equilibrium. Overall, the same anisotropy was observed and the results are in good agreement, even if a higher value for C 33 was obtained in our case and may be the consequence of the consideration of perfectly planar molecules.

Fourth order elastic tensor dependence on T and P

In Figure 3.2 is represented the dependence of the elastic coefficients with pressure. In Table 3.3, it can be observed that the elastic constants are not subjected to a noticeable dependence on temperature. On the other hand, one can notice in Table 3.4 that they undergo a significant increase with pressure, but with no change in anisotropy. The ratios between the C ij at 10 GPa and the C ij at 0 GPa are almost the same for all C ij . The values presented in Table 3

Large deformations in molecular dynamics

Prescribed deformation paths

Due to the high anisotropy of TATB, the mechanisms of irreversible deformation are suspected to strongly depend on the direction of loading.

The understanding of plastic activity and mechanical behavior of TATB single crystal is important at the atomic scale if one wants to build a constitutive law at larger time and space scales. To be able to detect all possible mechanisms, we propose a methodology to perform large deformation 3D-periodic MD simulations. In order to apply any type of deformation (pure shear, pure compression, compression with lateral distension, etc.) along any direction, we prescribe a velocity gradient

L = Ḟ • F -1
defined by a deformation gradient tensor F (t) and its time derivative Ḟ . The deformation gradient tensor F is constructed as follows:

F = I + α(t)m ⊗ m + β(t)n ⊗ n (3.18)
where α and β are time dependent functions, ⊗ is the dyadic product, (m, n) are the vectors used for the construction of F as described in For each MD simulations discussed in this thesis, a non-orthogonal basis (a, b, c) defines the crystal lattice. A second orthogonal basis (x, y, z) is defined by using the convention such that a is aligned with x, b is lying in the (x, y) plane and c is oriented in the upper cartesian half-space and makes angles of β and α with a and b, respectively. For MD simulations that involve time dependent prescribed deformations, a 100 ps trajectory was calculated in the NVT ensemble in order to reach equilibrium, using cell parameters obtained from a Parrinello-Rahman isostress-isothermal (N σT ) simulation (at the same temperature and pressure). The final configuration is then used as initial conditions for subsequent MD simulations with deformations (i.e. t = 0 for the time dependent deformation gradient tensor F (t)). Then, the MD supercell is dynamically strained, knowing its average deformation through the application of F to the (a, b, c) system's boundaries only (see Figure 3.4). A velocity field is applied at initial time before deformation (i.e. after 100 ps of equilibration)

and is consistent with the velocity gradient tensor L deduced from F (i.e.,

v(t = 0) = L • x(t = 0)).
No rescaling of the atomic positions are used in this method, which allows for a consistent motion of atoms/molecules regarding the applied strain and avoids the generation of acoustic waves in the material.

Local gradient deformation tensor computation

The aformentioned deformation, applied by dynamic modification of cell parameters, may lead to an irreversible evolution of the system through the onset of local defects. Two different analyses are available to capture these mechanisms.

When full dislocations (nucleation and motion) form, the original crystal structure is recovered behind the moving dislocation, but is strongly perturbed at the dislocation core. Local detection of loss of crystalline structure (as the one proposed, e.g., in the DXA algorithm [START_REF] Stukowski | Extracting dislocations and nondislocation crystal defects from atomistic simulation data[END_REF]) For these reasons, the computation of the local Lagrangian deformation gradient tensor F with respect to a reference defect-free state is preferred in the following. This measure has been proposed by several authors (Gullett et al. [START_REF] Gullett | A deformation gradient tensor and strain tensors for atomistic simulations[END_REF], Tucker et al. [START_REF] Tucker | Continuum metrics for deformation and microrotation from atomistic simulations: Application to grain boundaries[END_REF], Denoual et al. [START_REF] Denoual | Phase-field reaction-pathway kinetics of martensitic transformations in a model fe 3 Ni alloy[END_REF]) and more recently by Kober et al. [START_REF] Kober | Characterizing atomistic geometries and potential functions using strain functionals[END_REF] based on strain functionals.

highlights
In continuum mechanics, the current material points' position at time t in the strained state refers to the position in a reference state through the deformation gradient tensor F (t) with x(t) = F (t) • X, where

x and X are the current and initial positions, respectively. We approximate the deformation gradient tensor F (t) by a tensor F i defined for each molecule by a least square method that uses the relative displacement of its neighboring molecules (see appendix A.3 for further details).

Atomic positions filtering

Additionally, we used in our MD simulations with prescribed deformation paths an atomic position filtering based on the following first order differential equation:

r + t c ṙ = r (3.19)
where r and r are the filtered and real atomic positions respectively with t c a parameter with dimension of time. By considering ṙ = r(t+1)-r(t) ∆t , the resolution of this differential equation leads to the following expression for the filtered atomic positions:

r(t + 1) = αr(t) + (1 -α)r(t) (3.20)
where t is the time and α = ∆t tc with ∆t the MD simulation time step. The parameter t c is material-dependent and is computed from material properties using the Debye frequency to filter out phonon vibrations. For TATB, this leads to a value α = 0.01. This filtered position is used as an output only and has no influence on the modelled physics. This method will be used in chapter 5.

Peierls-Nabarro-Galerkin model for dislocations

The Peierls-Nabarro-Galerkin (PNG) method is a technique developed to calculate complex dislocation core structures, e.g., when the core is spread along different non-parallel glide plane [START_REF] Denoual | Modeling dislocation by coupling peierls-nabarro and element-free galerkin methods[END_REF], with complex crystallographic structure [START_REF] Metsue | Peierls-nabarro modelling of dislocations in diopside[END_REF], or submitted to a dynamic loading [START_REF] Denoual | Modeling dislocation by coupling peierls-nabarro and element-free galerkin methods[END_REF]. These calculations are usually performed at 0 K, which is consistent with the definition of the dislocation Peierls stress, the stress needed to overcome the lattice resistance at 0 K. One can notice however that the vibration frequency of a dislocation with a large core (as found in TATB, see chapter 4) is several orders of magnitude lower than the frequency of atomic vibrations. A time-scale separation could thus be introduced between the dislocation core structure and the molecules positions and orientations at equilibrium. This opens the way to the definition of a γ-surface (i.e. the generalized stack-ing fault energy E gsf (f ) that represents the displacement jump cost for a displacement f on a glide plane, see Figure 3.5) at finite temperature that can be used in a PNG dislocation core structure calculation. where E isf is the inelastic stacking fault energy (also called the "crystalline energy"), ε 0 is the stress-free strain and B is the Birch coefficients (or stress/strain) tensor [START_REF] Birch | Finite elastic strain of cubic crystals[END_REF] which is deduced from the elastic tensor C with the following formulae [START_REF] Wallace | Statistical Physics Of Crystals And Liquids[END_REF][START_REF] Wallace | Thermoelasticity of stressed materials and comparison of various elastic constants[END_REF] :

b 0 1 ||η|| x b L ε = η(x)ε 0 ε 0 = cste ε 0 = ||b|| 2L (n • b + b • n) PNG n x
B ijkl = C ijkl -P(δ jl δ ik + δ il δ jk -δ ij δ kl ), (3.22)
where P is the pressure (P > 0 for compression).

The determination of E isf depends on material properties and governs the inelastic behavior of the dislocation core. This energy is built such that the computed total energy from the continuum fields exactly equals E gsf (from MD), leading to a minimization relation connecting these two energies [START_REF] Denoual | Dynamic dislocation modeling by combining peierls nabarro and galerkin methods[END_REF]. Thus, the only input here are the generalized stacking fault energies E gsf obtained by independent MD calculations of block-sliding configurations (that will be detailed in chapter 4), and the complete fourth-order stiffness tensor C obtained by MD calculations of homogeneous deformations (presented in this chapter). An interpolation is performed on the γ-surfaces (inferred by MD simulations) in order to obtain a much more refined representation of E gsf . The interpolation error has been set to 0.5 mJ.m -2 , leading to a very accurate discretization of the γ-surface.

Conclusion

In the following, classical molecular dynamics simulations are performed with the Bedrov force field and rigid molecules to investigate the mechanical behavior of TATB single crystal under various conditions in temperature and pressure. In chapter 4, the dislocation core structures in TATB basal planes are calculated in temperature and pressure with the Peierls-Nabarro-Galerkin model introduced in section 3.3. Then, the mechanical response of TATB single crystal under directional loading is studied in chapter 5 using the method presented in section 3.2. Finally, the results of chapters 3, 4 and 5 are gathered in order to build the early bricks of the mesoscopic model and to study the mechanical behavior of both TATB single and polycrystal under triaxial compression in chapter 6.

Chapter Summary

• The force field used in this thesis has been validated through the computation of TATB cell lengths a, b, c, angles α, β and γ, and its fourth order elastic tensor C.

• TATB is found to be very stable with hydrostatic pressure, and both cell parameters and elastic constants are in good agreement with the results from the literature.

• The anisotropy, reflected by the C components is preserved when the material is subjected to hydrostatic pressure.

• The consideration of rigid molecules in our MD simulations is thought to be a reasonable hypothesis.

• The prescription of deformation along paths makes possible the exploration of a large range of non-hydrostatic pressures of any kind, i.e. uniaxial and triaxial compressions, pure shear etc.

• Such deformations may lead to the nucleation of local defects, the identification of which is necessary to understand the material behavior. Thus, the computation of the local Lagrangian deformation gradient tensor F was performed at each molecular COM position with respect to a reference defect-free state.

Résumé substantiel en franc ¸ais

Les méthodes de simulation utilisées durant la thèse sont présentées dans ce chapitre. La méthode de la dynamique moléculaire classique (DM) est présentée en partie 3.5.1 avec des étapes de validation du potentiel utilisé. Ensuite, en partie 3.5.2 nous présentons la méthode des chemins de déformation en DM. Enfin, une introduction au modèle de Peierls-Nabarro-Galerkin est proposée dans la partie 3.5.3.

d dt ∂L ∂ qi - ∂L ∂q i = 0, i = 1, ..., N, (3.23)
ce qui constitue une base pour la seconde loi de Newton. Généralement, la fonction Lagrangienne s'écrit comme suit :

L = K(q, q) -U (q), (3.24)
avec K et U des fonctions d'énergie cinétique et potentielle. On peut également raisonner en terme de Hamiltonien H avec :

H = i qi • p i -L = K + V, (3.25)
o ù l'impulsion généralisée p i est définie par p i = ∂L ∂ qi

. Cette formulation mène aux équations du mouvement de Hamilton :

qi = ∂H ∂p i , ṗi = - ∂H ∂q i . (3.26)
Ces équations ne sont en général pas résolvables analytiquement et sont discrétisées sur des intervalles de temps ∆t très courts de l'ordre de la femtoseconde. Lorsque les interactions entre particules sont régies par un potentiel inter-atomique U (q 1 , q 2 , ..., q i , ..., q N ), chaque particule est impactée par le mouvement de ses particules voisines au temps t et la force F i agissant sur la particule i est calculée comme :

F i = -∇U , (3.27)
o ù ∇ est l'opérateur gradient.

Molécules rigides

Certaines interactions intra-moléculaires (comme les effets de vibrations des liaisons) ont des temps caractéristiques (≈ 0, 1 fs) inférieurs aux interactions inter-moléculaires selon le pas de temps choisi pour les simulations de DM (≃ 1, 0 fs dans notre cas) et peuvent être exclues afin de considérer les molécules comme rigides. Les trajectoires des molécules sont alors calculées à partir de leur centre de masse (CDM) et de leur orientation. Le formalisme des quaternions introduit par Evans [START_REF] Evans | On the representatation of orientation space[END_REF] permet de modifier le Lagrangien afin de tenir compte de ces rotations (voir annexe A.1). Dans ce cas, le Lagrangien s'écrit :

L = 1 2 mV 2 + 1 2 Iω 2 -U , ( 3
Q NVE = 1 N !h 3N δ E(q, p) -E 0 dqdp, (3.30)
o ù h est la constante de Planck, dqdp une intégrale dans l'espace des phases, E(r, p) l'énergie de la configuration, E 0 l'énergie initiale et δ le symbole de Kronecker.

Ensemble canonique

L'ensemble canonique NVT n'est pas à énergie fixe mais à température T fixe. Le système échange donc de l'énergie avec un thermostat qui impose une température de consigne. La fonction de partition de cet ensemble thermodynamique est non uniforme et suit une distribution de Boltzmann :

Q NVT = 1 N !h 3N e -E(q,p) k b T dqdp, (3.31)
o ù k B est la constante de Boltzmann et E l'énergie du système. Les méthodes les plus connues permettant de coupler le système avec un thermostat sont celles de Berendsen et al. [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF], de Nosé-Hoover [START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF] et de Langevin [START_REF] Langevin | Sur la théorie du mouvement brownien[END_REF].

Dans ce travail, toutes les simulations de DM ont été effectuées avec le code STAMP [START_REF] Soulard | Molecular dynamics study of the micro-spallation[END_REF]. Les simulations ont été lancées sur le supercalculateur TERA1000 au CEA. Les représentations et différentes images ont été effectuées à l'aide de logiciels tels que Visual Molecular Dynamics Software, ParaView et Blender.

Potentiel pour le TATB avec molécules rigides

Bedrov et al. [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF] ont adapté un potentiel développé initialement par Borodin et al. [START_REF] Borodin | Polarizable and nonpolarizable force fields for alkyl nitrates[END_REF] pour modéliser le TATB. L'expression de ce champ de forces est la suivante :

U tot = U NB (r) + bonds U bond (r ij ) + bends U bend (θ ijk ) + dihedrals U dihedral (φ ijkl ) + improper dihedrals U imp (φ imp ijkl ) (3.32)
o ù U NB est une fonction régissant les interactions entre atomes non liés.

Ce potentiel consiste en une somme d'un terme de répulsion-dispersion à deux corps, d'interactions coulombiennes et d'une énergie de polarisation :

U NB (r) =U RD (r) + U coul (r) + U pol (r) = i>j A αβ e (-B αβ r ij ) -C αβ r 6 ij + D 12 B αβ r ij 12 + i>j q i q j 4πǫ 0 r ij - 1 2 i µ i E 0 i . (3.33)
o ù µ i est le dipôle induit au centre de force i, ǫ 0 la permittivité du vide et E 0 i le champ électrique généré par les charges partielles. A αβ et B αβ sont des paramètres de répulsion tandis que C αβ est un terme de dispersion entre atomes i et j de type α et β. Le terme D(12/(B αβr ij )) 12 , avec D = 5 × 10 -5 kcal mol -1 pour n'importe quelle interaction de paire, est proche de zéro pour des séparations atomiques mais devient prédominant pour r ij < 1 Å, assurant que U RD (r) est répulsif pour des courtes distances. Pour des interactions hétéro-atomiques, les règles de mélange de Waldman-Hagler [START_REF] Waldman | New combining rules for rare gas van der waals parameters[END_REF] sont utilisées :

A ij = A ii A jj B 6 ij B 3 ii B 3 jj ; B ij = 2 B ii -6 + B jj -6 1 6 ; C ij = C ii C jj . (3.

Discussion sur le choix des molécules rigides

Dans un travail récent [START_REF] Pal | Contribution of molecular flexibility to the elastic-plastic properties of molecular crystal α-rdx[END_REF], l'influence du degré de flexibilité des molécules du α-RDX sur ses propriétés élasto-plastiques a été étudié.

Les auteurs ont montré que les propriétés peuvent être très affectées par certaines hypothèses de rigidité de la molécule. Cependant, la molécule de RDX possède plusieurs configurations conformationnelles ce qui n'est pas le cas de celle de TATB. En effet, le noyau benzenique lui confère une grande planéité. C'est pourquoi nous avons choisi d'effectuer nos simulations de DM en utilisant des molécules rigides, ce qui permet un gain en temps de calcul considérable (10 fois plus rapide).

Équation d'état

Des simulations de Parrinello-Rahman [START_REF] Parrinello | Crystal structure and pair potentials: A molecular-dynamics study[END_REF][START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] ont permis d'obtenir l'évolution des paramètres de maille du TATB en température et en pression avec des constantes de couplage fixées à 1, 0 fs et 0, 1 ps, res-pectivement. Une boîte de simulation 3D-périodique de 5a × 5b × 7c avec les paramètres expérimentaux de Cady et Larson [START_REF] Cady | The crystal structure of 1,3,5triamino-2,4,6-trinitrobenzene[END_REF] contenant 350 molécules a été utilisée. Les paramètres de maille du TATB pour chaque jeu de (T , P ) ont été obtenus en moyennant leurs valeurs sur les dernières 150 ps d'un système équilibré et sont donnés dans le Tableau 3.2.

Les angles α, β et le paramètre c dévient de quelques pourcents des valeurs de Cady et Larson [START_REF] Cady | The crystal structure of 1,3,5triamino-2,4,6-trinitrobenzene[END_REF], ce qui est certainement d û à l'hypothèse de planéité des molécules. En Figure 3.1 est représentée l'évolution des paramètres de maille du TATB en fonction de la pression ainsi qu'une équation d'état. Ces résultats sont comparés avec différentes valeurs expérimentales [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF][START_REF] Olinger | The hydrostatic compression of explosives and detonation products to 10 gpa (100 kbars) and their calculated shock compression : Results for petn, tatb, co2 and h2o. 6[END_REF][START_REF] Plisson | Equation of state of 1,3,5-triamino-2,4,6trinitrobenzene up to 66 gpa[END_REF] ainsi qu'avec des données DFT [START_REF] Pineau | Cold curve of tatb[END_REF]. Nos résultats sont en bon accord avec la littérature bien qu'un léger écart soit perceptible à pression ambiante. Puisque cette thèse se focalise principalement sur la réponse mécanique du TATB, nous estimons que les légères différences structurelles en comparaison avec la littérature [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF][START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF][START_REF] Cady | The crystal structure of 1,3,5triamino-2,4,6-trinitrobenzene[END_REF] n'auront pas une grande influence sur les conclusions générales concernant le comportement du TATB. [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF][START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] bien qu'une différence notable soit présente sur C 33 qui est plus élevée dans notre cas, ceci étant certainement d û à la considération de molécules rigides parfaitement planes.

Propriétés élastiques

Le

Influence de la température et de la pression

Les valeurs des constantes élastiques en fonction de la pression sont données dans le Tableau 3.4. Dans le Tableau 3.3, on note que l'élasticité est peu dépendante de la température pour la gamme étudiée. Concernant l'évolution en pression, le ratio entre les C ij à 10 GPa et 0 GPa reste relativement constant, confirmant le maintien de la forte anisotropie du TATB.

Grandes transformations en dynamique moléculaire

Chemins de déformation imposés

Afin d'identifier de nouveaux mécanismes de déformation pour le monocristal de TATB, nous proposons une technique de simulation permettant d'appliquer une déformation quelconque à une boîte de simulation de dynamique moléculaire. La sollicitation est prescrite par un tenseur gradient de la transformation F (t) et un taux de déformation L = Ḟ • F -1 . Le tenseur gradient de la transformation est défini comme :

F = I + α(t)m ⊗ m + β(t)n ⊗ n (3.
E = 1 2 V ∇u -ηε 0 B ∇u -ηε 0 dr énergie élastique + S E isf [η] dr énergie inélastique (3.36)
o ù E isf est l'énergie inélastique de faute d'empilement, ε 0 est la déformation en contrainte relaxée et B est le tenseur d'élasticité de Birch [START_REF] Birch | Finite elastic strain of cubic crystals[END_REF] qui est déduit du tenseur d'élasticité C via la formule suivante [START_REF] Wallace | Statistical Physics Of Crystals And Liquids[END_REF][START_REF] Wallace | Thermoelasticity of stressed materials and comparison of various elastic constants[END_REF] :

B ijkl = C ijkl -P(δ jl δ ik + δ il δ jk -δ ij δ kl ), (3.37)

CHAPTER 4 DISLOCATIONS IN TATB (001) BASAL PLANES

Introduction

The dislocation core structures in the (001) basal plane of TATB are investigated in this chapter. The objective here is to determine if dislocations can be considered as discrete lines as in common metallic alloys, and thus using standard crystal plasticity constitutive laws for TATB single crystal. In This obviously indicates that the most favorable glide planes that have to be considered are the P1 and P2 ones.

The main mechanism that controls plasticity in crystalline materials is known to be the dislocations but their existence in TATB has never been observed. The unusual elastic anisotropy, as well as the very low shear energies along the gliding directions, impose to consider possible large dislocation cores, and thus huge (and expensive) simulation setup. To overcome this difficulty, a numerical method derived from the Peierls-Nabarro (PN) model [START_REF] Bulatov | Computer simulations of dislocations[END_REF] called Peierls-Nabarro Galerkin (PNG, see [START_REF] Denoual | Dynamic dislocation modeling by combining peierls nabarro and galerkin methods[END_REF][START_REF] Denoual | Modeling dislocation by coupling peierls-nabarro and element-free galerkin methods[END_REF]), which has been presented in section 3.3, is considered here.

In this method, the dislocation core structure emerges through the competition of a potential energy deduced from the γ-surface and the elastic Parts of this chapter previously appeared in: energy. It has been used for numerous studies dedicated to the simulation of complex dislocation core structures for metals [START_REF] Denoual | Modeling dislocation by coupling peierls-nabarro and element-free galerkin methods[END_REF] or earth mantle materials [START_REF] Metsue | Peierls-nabarro modelling of dislocations in diopside[END_REF][START_REF] Giacomazzi | Dislocation properties of coesite from an ab-initio parametrized interatomic potential[END_REF].

For the latter studies, the PNG computations were performed at 0 K to provide an estimate of the so-called Peierls-stress [START_REF] Peierls | The size of a dislocation[END_REF], defined as the ultimate shear stress without thermal activation. Noting the important dependency of thermodynamical and mechanical properties with T , and keeping our focus on dislocation existence and core structure, we propose in the following a method to obtain, through a PNG calculation, an estimation of dislocation core structure under non-zero temperature.

The understanding of plasticity mechanisms in TATB single crystal is essential to perform simulations of plastic deformations at nanometer scale as done by Koslowski et al. [START_REF] Cao | Rate-limited plastic deformation in nanocrystalline ni[END_REF][START_REF] Lei | The role of partial mediated slip during quasi-static deformation of 3d nanocrystalline metals[END_REF] who studied the role of partial dislocations and stacking fault ribbons in 3D nanocrystalline materials depending on the grain size, using a phase field theory of dislocation dynamics [START_REF] Koslowski | A phase-field theory of dislocation dynamics, strain hardening and hysteresis in ductile single crystals[END_REF]. The methodology presented hereby is dedicated to the study of dislocations in TATB and their thermal and pressure dependence.

Gliding properties of the P1 and P2 (001) planes

γ-surfaces computation

The γ-surface is a potential energy surface calculated by prescribing a gliding f along the dense plane of a perfect crystal. An additional relaxation is usually considered in the direction perpendicular to the gliding surface (see [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] for a γ-surface in TATB at 0 K). Due to this relaxation, and consistently with common practices that use γ-surface to calculate dislocation core structure without thermal activation, this potential energy is calculated at 0 K only.

Due to the sensitivity of the TATB single crystal on temperature, the dislocation structure calculated at 0 K could differ from the one obtained in temperature. Since the vibration frequency of a dislocation with a large core is several orders of magnitude lower than the frequency of molecules vibrations, a time-scale separation is introduced, allowing the calculation simulations at 300 K to obtain the lattice parameters in pressure, the same structure (i.e. triclinic) is obtained but with a reduction of the interplanar distance. Thus, as shown in Figure 4.3, the γ-surface in pressure preserves its shape but the energetic barriers are higher. At 10 GPa, the energy maximum is almost 4 times larger than at 0 GPa. The γ-surfaces for P 1 and P 2 planes exhibit the same energy barriers, due to the same interplanar distance involved. However, due to the centrosymmetry of TATB single crystal, the two γ-surfaces are in fact rotated relative to each other.

Balance between elastic and inelastic energies

The γ-surfaces do not experience a strong shape evolution under pressure, consistently with the high stability of the triclinic structure of TATB and the only effect is an increase of the energy barriers. This can be checked on the Minimum Energy Path (MEP) plots computed for a gliding of a [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] and b[010] by using the string method [START_REF] Weinan | String method for the study of rare events[END_REF] for different pressures (see Figure 4.4). In these MEP graphs, the energy values have been normalized by an estimation of the average shear modulus within We have considered planar and rigid molecules for our MD simulations in order to compute elastic constants and γ-surfaces. Results presented in the literature for both stiffness tensor [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF][START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] and γ-surfaces [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] have been obtained through MD simulations with flexible molecules.

Since our results with rigid molecules are similar to those obtained previously [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] with flexible molecules both for γ-surfaces and elasticity, we consider our hypothesis to be justified, thus allowing to use these results

as input of an upper-scale simulation such as a dislocation core structure computation with the PNG model.

Dislocation cores structure

For all calculations, a 2D mesh is considered with dimensions of 1024 TATB unit cells in the x direction (that corresponds to crystallographic direction a), 256 unit cells in the direction perpendicular to the glide plane (direction z), and one unit cell in the y direction. The nodal resolution is equal to two Element Galerkin nodes per lattice edge in x and y directions (i.e., 2048 × 512 nodes are considered). Last top and bottom nodes rows have a displacement constrained in the y direction, and have no constraint in the x direction. A discrete dislocation is introduced at the volume center and its structure relaxed to its minimal energy following the procedure proposed by Denoual [START_REF] Denoual | Modeling dislocation by coupling peierls-nabarro and element-free galerkin methods[END_REF]. The dislocation structure is then analyzed in terms of the spreading of the η field representing the inelastic part of the disregistry [START_REF] Denoual | Dynamic dislocation modeling by combining peierls nabarro and galerkin methods[END_REF].

Results for the dislocation core at 0 K and ambient temperature for edge and screw dislocations along [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] and [010] directions, for the two glide planes are first presented. Then, only the P 1 plane is discussed to illustrate the impact of pressure on edge and screw dislocations along a and b, since the behavior in the two glide planes P 1 and P 2 is very similar regarding their energy landscape and dissociation paths.

[100](001) dislocations

The dislocation core structure is determined by analyzing the Frobe- First, it appears that the [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] edge dislocation gliding on the (001) plane leads to a very complex core structure which spreads out over around 125 Burgers vectors at 0 K. The [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] screw dislocation also exhibits a complex (but thinner) core, spread over 75 Burgers vectors. Not only those widths are several orders of magnitude higher than the one that could lead to a measurable dislocation friction, but they also reveal that dislocations in TATB cannot be considered as discrete lines but rather as stacking fault ribbons.

nius
The second set of simulations shows that an increase of T leads to a reduction of the spreading length for both edge and screw dislocations as shown in spreading within the glide plane is then expected. Here as well, screw dislocations spread over a smaller distance than edge ones.

The results for the P 2 plane depict similar trends. For both edge and screw dislocations, for the two glide directions and the two glide planes, dislocations under ambient pressure have no localized core, and setting the temperature at 300 K sligthly reduces the core spreading about approximatively 25 %.

Effect of pressure on the dislocation cores

Because the PNG model consists in finding the equilibrium between elastic and stacking fault energies, one can anticipate how the dislocation core structure is affected by the pressure. All dislocations seem to split into several partials.

In order to analyze the evolution of dislocation core with pressure, we defined ∆ SL as the distance between the two most distant peaks of |∇η(x)| which represents the equilibrium distance between partial dislocations.

Conclusion

An original method to compute γ-surfaces and elastic constants in temperature and pressure through MD simulations has allowed the investigation of dislocation core structures in TATB molecular crystal.

Elastic constants dependence in pressure shows that anisotropy is not enhanced as the ratio between C ij remains constant as pressure increases.

The γ-surfaces energy landscapes are sligthly smoothed with temperature. Their shape is preserved between ambient pressure and 10 GPa

Chapter Summary

• A new method for the computation of γ-surfaces in temperature has been introduced. The γ-surfaces of TATB (001) planes exhibit very low energies, indicating that dislocations in these planes could easily glide.

• The energies along the minimum energy paths on the γsurfaces for full Burgers vectors, normalized by an estimation of the shear modulus in the (001) plane, show that the dislocation core structure is maintained with pressure.

• Both elastic constants and γ-surfaces have been used as an input in a Peierls-Nabarro-Galerkin model, in order to compute the first dislocation core structures in TATB (001) planes.

• We show that [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF](001) and [010](001) edge and screw dislocations split in three or more partials. For edge components, the distance between the two most distant partial dislocations ∆ SL is 1.4 times greater than for screw components.

• Dislocations in TATB basal planes can not be considered as discrete lines, but rather as large stacking fault ribbons. Dislocation friction on the lattice in these planes in unlikely to play a role for the viscoplastic behavior of TATB material.

Résumé substantiel en franc ¸ais

Les structures de coeur de dislocations dans le plan basal (001) du TATB sont étudiées dans ce chapitre. Une vue dans la direction [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] du monocristal de TATB est proposée en

Nucleation stress surface for TATB single crystal

Simulation characteristics

MD simulations were performed with a 3D-periodic 10a × 10b × 14c simulation cell containing 2800 molecules. Cell parameters at different temperatures and pressures were previously obtained through Parrinello-Rahman simulations. At first, the system is equilibrated in the NVT ensemble. Then, the large strain methodology presented in section 3.2.1 is used to apply directional pure shear to the simulation box. The plane defined by the two vectors m and n is chosen to be perpendicular (see Figure 3.4) to the basal plane of TATB. Indeed, due to the large steric hindrance present in the molecular layers, we estimate that a pure shear with m and n lying in the layers is not likely to induce defects in these directions such as dislocations for example. Furthermore, with the choice made on the orientation of m and n, pure shear loadings are applied in directions which are expected to follow the TATB single crystal anisotropy.

Two conditions are applied concurrently for these simulations: volumepreserving strains are applied at constant strain rate, in order to construct the following expression for the deformation gradient tensor, function of time t: This mechanism is investigated and explained in detail in section 5.3.

F = I + (e At
The zone related to the lowest value of σ v (label 2) corresponds to the deformation paths where φ = 45 • . The deformation is equivalent to a pure shear (see Equation 5.1) in the basal plane (001). Glide is activated for stress below computable accuracy, owing to the weak van der Waals interactions between layers, which leads to the low energy barriers observed in the γ-surfaces [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF][START_REF] Lafourcade | Dislocation core structure at finite temperature inferred by molecular dynamics simulations for 1,3,5-triamino-2,4,6-trinitrobenzene single crystal[END_REF]. Due to the low friction in the TATB basal plane, no dislocations are nucleated, which is consistent with the main conclusion of chapter 4 on dislocations core structures for these directions.

Finally, in the zone 3 that looks like a "donut", the deformation is mainly compressive in the basal plane, inducing an alternating disorientation of molecular layers relative to their initial planar configuration (see Figure 5.6). This mechanism has already been observed in oriented block copolymers [START_REF] Makke | Nanoscale buckling in lamellar block copolymers: A molecular dynamics simulation approach[END_REF] and is referred to as "buckling". An analysis of this phenomenon is provided in section 5.4. It is important to note that non-basal plane dislocations and buckling deformation mechanism are triggered for similar values of the nucleation von Mises stress at standard ambient conditions.

Non-basal dislocations

As introduced in section 5.2, when pure shear is applied to the TATB single crystal with φ > 45 • , gliding of the basal plane is not activated but the overall deformation is mainly accommodated by dislocations nucleation in non-basal planes. The only way dislocations have to glide in transverse planes is through a hydrogen-bond breaking phenomenon, since inter-molecular bonding in the layers is governed by hydrogen bonds.

An example of dislocation propagating in the (0 11) plane is provided in The two slip systems in the (001) planes have already been studied, first by Mathew and Sewell [START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF] who computed the corresponding γ-surfaces and in this work [START_REF] Lafourcade | Dislocation core structure at finite temperature inferred by molecular dynamics simulations for 1,3,5-triamino-2,4,6-trinitrobenzene single crystal[END_REF] where the first estimate of the dislocation core structure was reported, as described in chapter 4.

The planes (100) and (010) are of similar nature and can be studied together. Among all possible Burgers vectors for dislocation nucleation within (100) and (010) planes, [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF], [010] and [001] Burgers vectors were never observed and we think that this type of gliding can be dismissed, due to the large steric hindrance present within the molecular layers. However, the gliding with a [001] component (thus implying hydrogen-bond breaking) in transverse planes has been observed in our simulations and occurs via formation of a complex stacking fault structure that entangles the planes P 1 and P 2. This very particular structure leads to a noticeable dilatancy within the dislocation core up to 10 -15 %, as explained in the next section.

[1 0 1] [1 0 1] a c b β [0 1 1] [0 1

Local dilatancy process

In undeformed TATB crystal, the packing of molecules is hexagonal

Activated slip systems

The procedure used to prescribe deformation paths in MD simulations, described in section 5. tensor. All the slip systems that have been found are listed in Table 5.1

(slip plane normal + Burgers vector).

Twinning-Buckling mechanism

For deformation paths that contain a compressive component in the basal plane, a chevron-like structure appears as shown in Figure 5.6.

Such a deformation can be seen as a composition of a pure rotation and a pure shear, which, in the case of TATB single crystal, systematically involves a homogeneous gliding in the (001) plane. It is not clear however if this mechanism is twinning (thus involving a well defined deformation that inverses the lattice) or buckling, which involves only a gliding slip system slip plane normal Burgers vector to characterize the behavior under compression of TATB single crystal in terms of stability, size and strain rate dependence. For the sake of simplicity, in the following we call this mechanism "Twinning-Buckling".

[ 1 0 1]( 1 0 1) ( 1 0 1) [ 1 0 1] [1 0 1](1 0 1) (1 0 1) [1 0 1] [0 1 1](0 1 1) (0 1 1) [0 1 1] [0 1 1](0 1 1) (0 1

Perfect twinning

The triclinic TATB molecular crystal is centrosymmetric with space group P1 [START_REF] Cady | The crystal structure of 1,3,5triamino-2,4,6-trinitrobenzene[END_REF], and consists of alternating stacking (ABABAB) of molecular layers. The hexagonal-like basal plane symmetry and stacking of the layers lead to a perfect twinning structure that can be obtained (at constant

Homogeneous shear deformation

The energetics along the twinning path is investigated in this section, with a particular focus on possible metastable phases (or twin structures)

and the stress needed to nucleate these phases. The twinning deformation is a simple shear in the basal plane (001), so that a comprehensive calculation of this energetics can be done by calculating the excess energy after an instantaneous deformation given by:

F = I + 1 d int f ⊗ n (5.2)
where f is the shear vector defined as

f = f 1 a[1 0 0] + f 2 b[0 1 0] with f 1
and f 2 varying from -0.5 to 0.5. n = z is the normal vector to the basal plane, d int the interplanar distance between molecular layers, assumed constant here, and I is the identity tensor.

This landscape must be calculated at 300 K to be consistent with the MD simulations done in this study. However, prescribing an overall deformation is not compatible with the scheme to integrate the equations of motion since an evaluation of the energy for each fraction of the shear vector f is required, without energy minimization. Assuming ergodicity for the system, the system is sampled over the equilibrium state by considering a large simulation supercell of 20a × 20b × 28c equilibrated in the NVT ensemble during 100 ps (with cell parameters obtained with the Parrinello-Rahman simulations at 300 K). Then, an instantaneous deformation through F is applied in order to evaluate the energetics along the prescribed deformation. A similar procedure has been proposed in chapter 4 to estimate the TATB basal γ-surface at 300 K [START_REF] Lafourcade | Dislocation core structure at finite temperature inferred by molecular dynamics simulations for 1,3,5-triamino-2,4,6-trinitrobenzene single crystal[END_REF]. and twinned configurations such as the force derived from the potential is tangent to the MEP. They are calculated with the string method [START_REF] Weinan | String method for the study of rare events[END_REF].

Behavior under pure compression

The deformations prescribed in this section are simple-compressive in the basal plane, and therefore different than the isochoric deformation used in section 5.2, (where it was applied to the study of the onset of plasticity under shear loading). Thus, we apply pure compression to the molecular layers with φ = 0 • and the deformation gradient tensor is constructed from the vector m only.

An orthonormal basis (x, y, z) is defined with x parallel to the triggered mechanisms.

Two types of simulations were undertaken: the first consisted in applying an instantaneous compression ratio R followed by an equilibration of the system energy in the NVT ensemble, devoted to observe the arrangement of the molecular layers and its dependence on the compression ratio R and direction m. The second type of simulation involved a first step of equilibration in the NVT ensemble for 50 ps followed by a constant strainrate deformation up to a chosen compression ratio R = 5% which permits to study the strain-rate dependence of the generated microstructure.

Effect of strain rate

The case of prescribed deformation is considered here for several constant strain rates, corresponding to a final strain of 5% with trajectory lengths between 50 ps and 500 ps. In ǫ, E being a strain measure with E = 1 2 (F T F -I) the Green-Lagrange strain tensor computed from the local deformation gradient tensor for each molecule. The increase of strain rate leads to higher nucleation stresses before stress relaxation due to a large number of nucleation sites, compared to what is observed at lower strain rates. This is consistent with a hypothesis of a nucleation rate driven by thermal activation, with the activation energy decreasing with increasing local stress. Each nucleation site produces a chevron that will relax the stress in a zone increasing with time. Such competition between nucleation and screening is usually well described by the Avrami model [START_REF] Avrami | Kinetics of phase change. i general theory[END_REF] that predicts a decrease of the characteristic length for an increasing nucleation rate, as observed in Figure 5.11.

In the same way, the size of the large deformation bands (denoted in the following as "chevrons") in the microstructure is also related to the strain rate, with larger chevrons at lower strain rates. For the lowest strain rate, we observe a very progressive growth of few chevrons toward a typical converged length of about 50 -60 TATB lattices. However, this length is certainly a finite-size effect since when the entire interface between chevrons becomes perpendicular to the loading direction, the remaining chevrons are disconnected and their merging stops. Thus larger system size would lead to larger chevrons.

Dependence on loading direction

Multiple directions of compression were investigated for this study at the same strain rate of 2 × 10 8 s -1 . The TATB supercell was compressed along the direction m = (cos θ, sin θ, 0) for θ = 30, 75 and 90 • .

Chevrons are nucleated with interfaces perpendicular to the compression direction as shown in Figure 5.12. Snapshots of three molecular dynamics simulations (different rows) are represented, corresponding to different compression angle θ (30, 75 and 90 • ), and intermediate strains (different columns) ǫ (3, 4 and 5%). Starting from an equilibrated and defect-free single crystal, the trigger of the first disorientations happens around 3% strain, with a homogeneous nucleation of tiny chevrons. As the deformation progresses, a hierarchical evolution is observed with the merging of chevrons to larger structures. In this case, the TATB sample size was 100 × 100 × 10 lattices. Once the structures were equilibrated, interfaces were constrained by the 3D periodic cell and the growth could not continue.

Microstructure evolution

The buckling involves significant shear of the basal plane. In this section, the role of the energy landscape calculated for simple shear in sections 5.4.1 and 5.4.2 on the formation of the microstructure is assessed.

As a first step, an estimation of the buckling threshold as the highest defor- 92% molecules of the system, i.e. more than 92% of the system is affected by a precise and well-defined deformation.

In Figure 5.13, the simplified histogram calculated at the end of the MD simulations has been reported, for compressions with θ = 0 , 30 , 75 and 90 • and a final strain equal to 5%. The histograms sets are represented by filled circles whose area is proportional to the fraction of the molecules in this state. For direction θ = 0 • , the deformation is only in C 13 and the sests are aligned with the direction of loading and centered on the strain-free minimum. In a similar way, when θ = 90 • , the two sets are also aligned with the direction of loading but most of the deformation is in C 23 . Concerning the case when θ = 75 • , the furthest set from the strain-free minimum is splitted into two subsets with one falling near the set of direction with 90 • . These observations show that when the direction of compression is aligned with the two minima (i.e. θ = 90 • ), the microstructure is determinated by twinning, with a strong influence of the underlying energy landscape. Conversely, the more the direction of compression is aligned perpendicularly to the two minima (i.e. θ = 0 • ), the more the microstructure is dominated by pure elasticity. For the range of loading directions investigated in this study, we can say that for θ ∈ [75 • , 105 • ] a splitting occurs and twinning is coupled with buckling. In Figure 5.14 are reported the energies along straight lines passing by the minimum of the α-TATB phase on the twinning energy landscape. We directly see that for θ = 0 • (i.e. perpendicularly to the two minima) the elastic energy well (yellow curve) is quadratic in its vicinity. Conversely, for θ = 90 • (i.e. aligned with the two minima), two energy wells are present (black curve), giving the possibility to fall in the vicinity of the twinning minimum. Thus, the proportion of the system that is whether in pure buckling or in pure twinning is determined by the loading direction. We termed this "in-between" mechanism Twinning-Buckling in the following.

The sets fall in the vicinity of the deepest energy well, confirming the role of the underlying energy landscape on the structure. This can be understood by comparing the excess energy for the minimum energy paths (MEPs) between two minima (see section (5.4.2)), estimated at 25.4 MJ m -3 , and the elastic energy available prior to the onset of buckling. This latter quantity is the elastic energy to be released by the buckling by considering the elastic energy stored by a deformation of 2.75% (the buckling threshold T B ) applied under uniaxial compression:

E B = 1 2 C 11 T 2 B = 18
.9 MJ m -3 . The two energies are therefore similar, demonstrating that even a small compression of a few percents is sufficient to overcome the energy barrier leading to the twinned configuration.

Furthermore, if only one minimum had been present on the energy landscape, with a classical quadratic elastic well, two sets with deformations of opposite signs and exactly same amplitude would have been observed which is not the case here. This implies the necessity to consider the twinning energy landscape during the construction of a constitutive law that would aim to model the Twinning-Buckling mechanism. on the loading direction. Indeed, the significant increase of the nucleation stress for the onset of dislocations (zone 1 in Figure 5.2) with respect to the pressure is consistent with the dilatancy of the dislocation core that could hinder the nucleation process as the pressure increases. On the contrary, the buckling mechanism (zone 3 in Figure 5.2) undergoes a weaker evolution for its nucleation stress with pressure. Thus, in the presence of large hydrostatic pressure, the strength of the hydrogen bonds network as well as the increase of the steric hindrance would prevent the nucleation of transverse dislocations and buckling would be the prevailing mechanism.

Size and strain rate effects

In the simulation described above, system size effects may have impacted the nucleation stress threshold. This must be understood in order to be able to build a constitutive law at larger time and space scales.

The applied strain rate was very high compared to standard experiments conditions at an engineering scale but more appropriate compared to the local material strain rate in a shock rise. Additional MD trajectories were carried out at a lower strain rate in order to investigate the effect of strain rate on the threshold of the two main mechanisms described above. First, three trajectories have been performed while increasing the system's size. Taking the cell parameters at 300 K and 0 GPa, MD simula-tions were performed with a 3D-periodic 40a × 40b × 60c simulation cell containing 179200 molecules. Deformation paths were performed with θ = 90 • and for three different values of φ: 15 • , 75 • and 90 • . The strain rate was the same as for the simulations at (300 K, 0 GPa) (i.e. 10 8 s -1 ). The aim of these three simulations was to quantify the size effects on the nucleation stresses. The different values of the von Mises stress at the nucleation of the first defect are reported in Table 5.2. As anticipated, increasing the size of the system decreases the nucleation stress threshold.

10a × 10b × 14c 40a × 40b × 56c ratio φ( • ) σ v (GPa) σ v (GPa) (%) 15
10 8 s -1 10 7 s -1 ratio φ( • ) σ v (GPa) σ v (GPa) (%) 0 1.
By multiplying the system size by 4 in each direction, the von Mises stress at nucleation has been lowered by approximately 20%, regardless of the loading direction. Larger system sizes have been investigated with no significant effect on the nucleation stress.

In order to study the effect of the strain rate, additional trajectories were undertaken with the same size as before (i.e. 10a × 10b × 14c). Deformation paths were chosen with θ = 0 • and φ every 15 • in [0, π/2], as for the computation of the nucleation stresses surface in Figure 5.2. The results are shown in Table 5.3, where one can see that the lower strain rate reduces the nucleation stress threshold by 14%.

Buckling: A linear elasticity problem?

The buckling mechanism seems to be unavoidable for TATB when the crystal is subjected to compression in the direction of its basal plane. This is mainly due to its layered structure and contrasts between strong intra-and weak inter-layer molecular interactions, confering very anisotropic elastic and plastic properties, as it can be seen through its fourth-order elastic tensor C. In this section, we propose a first modelling of the buckling by using purely elastic response.

Even if in the triclinic case the stiffness tensor is supposed to have 21 independent components, an approximation is made in the following to consider TATB as an orthotropic material since its elastic tensor contains multiple components near 0. For an ideal orthotropic material, C is given by:

C =              C 11 C 12 C 13 0 0 0 C 12 C 22 C 23 0 0 0 C 13 C 23 C 33 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 55 0 0 0 0 0 0 C 66              (5.3)
Taking the TATB complete fourth-order tensor computed through MD simulations (see chapter 3) at 300 K and 0 GPa, the simplified or-thotropic C reads:

C ≃              49
             , ( 5.4)
The understanding of the mechanical response of such anisotropic materials is far from being complete and many studies are dedicated to it.

Especially, anisotropic layered nanostructures exhibit complex mechanical response when they are stretched perpendicularly or compressed in the direction to the layers: in block copolymers for example [START_REF] Makke | Nanoscale buckling in lamellar block copolymers: A molecular dynamics simulation approach[END_REF], buckling of the hard lamellae is observed leading to a "chevron-like" microstructure.

This buckling instability can be studied analytically using elasticity theory as proposed by Read et al. [START_REF] Read | The chevron folding instability in thermoplastic elastomers and other layered materials[END_REF] and more recently by Makke et al. [START_REF] Makke | Nanoscale buckling in lamellar block copolymers: A molecular dynamics simulation approach[END_REF].

The volume average elastic energy density of a strained homogeneous material is given by:

e = 1 2 C ijkl ε ij ε kl (5.5)
with ε ij the macroscopic strain and C ijkl the elastic components. In the following, the 2D case is considered where a sample is compressed along the y direction and stretched in the z direction where y and z are parallel and perpendicular to the TATB layers respectively. Using the Voigt notation, Equation 5.5 reads:

2 e = C 22 ε 2 22 + C 33 ε 2 33 + 2C 23 ε 22 ε 33 + 2C 55 ε 2 23 .
(5.6)

Read et al. [START_REF] Read | The chevron folding instability in thermoplastic elastomers and other layered materials[END_REF] proposed an original approach consisting in the description of buckling as a local stretch deformation followed by a rotation, as illustrated in

                   ∆ = (∇ y v z ) 2 + (∇ y v y ) 2 1/2 sin θ = ∇v z /∆ ε 22 = ∆ -1 ε 33 = (1 + ∇ z v z )(1 + ∇ y v y ) -∇ z v y ∇ y v z -1 ε 23 = ∇ z v y (1 + ∇ y v y ) + ∇ y v z (1 + ∇ z v z ) /∆ (5.8)
A small perturbation u imposing a displacement along z only is introduced. Therefore, the overall displacement field v is a function of the macroscopic deformation (ε yy and εzz ) and u components:

       v x = 0 v y = ε yy v z = ε zz + u z (5.9)
with u z = U 0 sin(ky) a sinusoidal perturbation of wave vector k. Using the expression of v into the System 5.8 allows then to obtain the local deformations ε 22 , ε 23 and ε 33 as functions of the macroscopic deformation εyy and εzz : MD simulations and concerns the pure shear in the basal plane where the energies are very low due to the weak van der Waals interplanar interactions. Thus, the measure of such sensible elastic constants may have lead to slight errors. To lower the theoretical threshold for buckling to the observed value of ≃ 3 %, C 55 needs to be reduced from 1.9 GPa to 1.4 GPa,

             ε 2 22 = ε2 yy + U 2 0 k 2 2(1+εyy) 2 εyy (1 + εyy ) ε 2 33 = ε2 zz - U 2 0 k 2 2(1+εyy) 2 εzz (1 + εzz ) ε 22 ε 33 = εyy εzz + U 2 0 k 2 4(1+εyy) 2 (ε zz -εyy ) ε 2 23 = U 2 0 k 2 2(1+εyy) 2 (1 + εzz ) 2
showing the sensitivity of such an analytical estimation and the caution with which such a result must be understood. Yet, this difference on C 55

for the elastic onset of buckling is within the computational error bars of the stress. Buckling instability nucleation process may thus be described by elasticity. However, the obtained non-symmetric microstructures have to be described by the complete energy landscape that characterizes the coupling between twinning and buckling.

Conclusion

The computation of the yield stress under directional loading for TATB single crystal through MD simulations allowed for the identification of three distinct behaviors of the material. Gliding in the basal plane under pure shear is very easy to activate and is related to the very low energies on the γ-surface in these directions. This is consistent with a previous study of dislocation core structures for this plane presented in chapter 4 and in a previous paper [START_REF] Lafourcade | Dislocation core structure at finite temperature inferred by molecular dynamics simulations for 1,3,5-triamino-2,4,6-trinitrobenzene single crystal[END_REF]. Non-basal plane dislocations nucleation is a new plasticity mechanism for the TATB single crystal. It is shown to involve hydrogen-bond breaking and complex dislocation core structures with modification of the hexagonal in-plane molecular packing (for partial dislocations), and a noticeable core dilatancy. At the stage of this work, the informations regarding the dislocation-mediated plasticity for TATB single crystal still remains poor. Even if enough slip systems have been identified to be used in a crystal plasticity framework at upper scales, a lot of data is still missing concerning TATB crystal plasticity. Indeed, in order to be used in a mesoscopic modeling, additional informations would be needed such as Peierls stress values for activation, interactions properties within different dislocations, velocity rules etc. Unfortunately, due to the triclinic cell of TATB, computations of such objects are not possible at the moment and will be discussed as a perspective in chapter 7. Finally, when the strain consists mainly in a compression of the molecular layers, a buckling phenomena is inescapably involved, leading to a modification of the microstructure in a succession of chevrons of opposite deformations, with interfaces always normal to the compression direction.

The calculation of a homogeneous twinning energy landscape has also been performed. It revealed an energy minimum corresponding to a perfect twinning of the TATB single crystal. The transition between buckling and twinning is shown to be dependent on this energy landscape and the direction of loading. Furthermore, the buckling is a phenomenom that is triggered by an elastic instability. This has been evidenced by a theoretical study allowing to analyze this inevitable feature for layered materials with strongly anisotropic elastic properties. In Figure 2.1 in chapter 2 was shown an optical microscopy micrograph of a TATB-based PBX microstructure. Large deformation bands were clearly visible with random orientations and may be twinning-buckling mechanism signatures. Even if this twinning-buckling mechanism results in an equilibrated microstructure with successive positive and negative deformations, it is quite similar to the mechanism that leads to kink bands (as observed in polycrystalline ice [START_REF] Piazolo | Effect of local stress heterogeneities on dislocation fields: Examples from transient creep in polycrystalline ice[END_REF] for example). TATB single crystal behavior under large strains is very dependent on the loading direction: out-of-plane loading leads to dislocations nucleation, in-plane loadings systematically produce twinning-buckling whereas basal plane pure shear leads to simple homogeneous gliding. All these elements are the necessary ingredients to build a representative constitutive law for the TATB single crystal behavior under large strains, which ideally would take into account dislocations nucleation, perfect lattice twinning and buckling mechanism.

Chapter Summary

• Prescribed deformation path MD simulations were performed to apply directional pure shear. The computation of the yield stress allowed for the construction of a 3D yield surface for TATB leading to three distinct behaviors.

• Pure shear in the basal plane provokes a gliding due to the low energies in these directions.

• MD simulations predicted nucleation of non-basal dislocations when the pure shear is applied in direction far from the basal plane. Dislocations involve hydrogen-bond breaking and complex core with a local dilatancy, that would hinder their nucleation under pressure. Four different glide planes were identified.

• When the strain involves a compression in the direction of the molecular layers, a buckling phenomenon is systematically observed, leading to a microstructure that consists in a succession of chevrons of opposite deformations, with interfaces normal to the direction of loading.

• A Twinning energy landscape has been calculated through the application of pure shear in the basal plane leading to a second energy minimum corresponding to a perfect twin. Buckling is shown to strongly depends on this energy landscape.

Résumé substantiel en franc ¸ais

Nous présentons dans ce chapitre une étude complète des mécanismes de déformation du TATB. La méthode présentée en partie 3.

Dislocations non-basales

En analysant la structure cristalline du TATB, on peut intuiter les différents plans dans lesquels des dislocations pourraient se déplacer. La

Mécanisme de maclage-flambage

Pour des déformations comprenant une compression dans les feuillets de molécules, une structure en chevrons apparaît (voir

T B = 2.75% est égale à E B = 1 2 C 11 T 2 B = 18, 9 MJ m -3 .
Les deux énergies sont assez proches, ce qui démontre que même une compression de quelques pourcents peut déclencher une microstructure comprenant des macles. De plus, si un seul minimum avait été présent dans le paysage énergétique, avec un puit d'énergie élastique quadratique, les couples de la déformation auraient été parfaitement centrés sur le minimum stable pour n'importe quelle direction, ce qui n'est pas le cas ici.

Effets de la pression, de la taille et de la vitesse de déformation

L'influence de la pression sur la surface de nucléation a été étudiée pour des pressions de (2, 5, 10 and 20 GPa). Afin de réduire le volume de calcul, nous avons choisi des trajectoires de déformation en cisaillement pur pour θ = 90 • et φ tous les 15 • dans [0, π/2]. Les résultats sont présentés en Figure 5.2. On voit que l'évolution de la contrainte de nucléation dépend du type de défaut activé et donc de la direction de sollicitation. En effet, la forte augmentation de la contrainte de nucléation pour les dislocations transverses est cohérente avec le phénomène de dilatance locale observé au coeur des dislocations, qui est limité lorsque l'on augmente la pression.

En revanche, le mécanisme de flambage est sujet à une évolution moins prononcée de sa contrainte de nucléation. Ceci pourrait indiquer qu'en présence d'une forte pression hydrostatique, le mécanisme de flambage serait favorisé.

Lors des simulations concernant la surface de contraine de nucléation, des effets de taille ont pu impacter la valeur du seuil. En augmentant la taille du système par 4 dans toutes les directions, la contrainte de nucléation, peu importe la direction de sollicitation, a été diminuée de 20%. De manière similaire, en diminuant la vitesse de déformation par 10, la contrainte de nucléation a été réduite de 14%.

Flambage : Une instabilité élastique ?

Une première étude analytique du flambage a également été proposée. En effet, cette instabilité est rencontrée dans les matériaux très fortement anisotropes soumis à des compressions uniaxiales dans les directions les plus rigides. Basée sur les travaux de Read et al. [START_REF] Read | The chevron folding instability in thermoplastic elastomers and other layered materials[END_REF], cette étude théorique mène à l'estimation d'un seuil en déformation à partir duquel le flambage peut être déclenché. Une fonctionnelle d'énergie élastique est obtenue en fonction des paramètres du matériau et des composantes de la déformation macroscopique qu'il subit (compression dans la direction rigide, traction dans la direction transverse). Cette énergie s'avère être très proche de zéro pour la configuration non déformée et négative dès que l'on comprime le matériau. Ainsi, pour une compression pure, le seuil théorique calculé pour le TATB est de l'ordre de 4%, representé par une étoile dans la Figure 5.17. De plus, une traction supplémentaire dans la direction transverse aura tendance à diminuer ce seuil de flambage.

Le seuil théorique de 4% est légèrement plus élevé que celui mesuré en DM. Ceci est peut-être d û aux simplifications effectuées sur le tenseur d'élasticité du TATB pour ce calcul analytique, prévu pour les matériaux de type orthotrope. En effet, les constantes de couplages proches de zéro dans le tenseur d'élasticité du TATB pourraient avoir un effet sur ce seuil théorique. Cette étude analytique du flambage nous a permis de montrer que ce mécanisme est principalement provoqué par une instabilité élastique. Cependant, la microstructure obtenue lors de nos simulations de DM est influencée par le paysage énergétique, qui conduit à une microstructure combinant maclage et flambage.

Conclusion

Dans ce chapitre, le calcul de la contrainte de nucléation pour un cisaillement pur en fonction de la direction de sollicitation a été effectué via des simulations de DM, menant à l'identification de trois comporte- non-basal dislocations are likely to nucleate for low pressures within four different identified slip systems, a perfect twinning of the TATB triclinic lattice is possible through a simple shear in the basal plane and small uniaxial compressions of the molecular layers systematically involves a buckling mechanism, with a resulting microstructure that strongly depends on the twinning energy landscape. The aim of this chapter is to present some preliminary results on simulations performed at the mesoscopic scale (µm), focusing on buckling.

When it comes to the construction of predictive constitutive law for shock loadings, the physical properties have to be evaluated within a certain range of pressure. In this chapter, we are particularly interested in Equations Of State (EOS) [START_REF] Birch | Finite elastic strain of cubic crystals[END_REF] that link volume and pressure. Such EOS are generally measured experimentally or calculated via numerical simulations (see [START_REF] Godwal | Equation of state theories of condensed matter up to about 10 tpa[END_REF] for a review) and then fitted with empirical functions (e.g. Birch-Murnaghan EOS [START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF][START_REF] Murnaghan | The compressibility of media under extreme pressures[END_REF][START_REF] Birch | Finite elastic strain of cubic crystals[END_REF], Vinet EOS [START_REF] Vinet | A universal equation of state for solids[END_REF]).

For very anisotropic materials, taking TATB as an example, the response to hydrostatic versus non-hydrostatic pressure leads to very different behaviors. Indeed, when subjected to hydrostatic pressure, these materials undergo different deformations according to their directional strength. Conversely, under triaxial compression, the deformation is directionally proportional, resulting in a non-hydrostatic stress. In opposition, for isotropic materials, these two types of loading lead to the same me-137 chanical response and no distinction is needed.

The chapter is organized as follows. In section 6.2 we introduce the two definitions of EOS that we consider in the following, that is when hydrostatic pressure is guaranteed on the sample (denoted HP EOS for Hydrostatic Pressure EOS), and when the material is subjected to a spherical deformation (denoted SD EOS for Spherical Deformation EOS). By spherical deformation, we mean a "homothetic" triaxial compression. Then, these two EOS are computed through molecular dynamics (MD) simulations and presented in section 6.3. In section 6.5, we first use an element free Galerkin method (EFG) to perform an elastic-buckling simulation (EB simulation) of the TATB single crystal under spherical deformation.

Finally, an EB simulation of TATB polycrystal is proposed and discussed.

Hydrostatic pressure vs. spherical deformation EOS for very anisotropic materials

The response of materials subjected to large deformations (e.g. during shock compressions) that involve the exploration of a wide range of pressures is of interest in many areas. A number of work has been undertaken to characterize the densification of matter with pressure, giving rise to various EOS e.g. the Vinet's EOS [START_REF] Vinet | A universal equation of state for solids[END_REF]:

P = 3K 0 x 2 (1 -x) exp 3 2 (K ′ 0 -1)(1 -x) , (6.1)
with x = V V 0 1/3 and K ′ 0 the pressure derivative of the bulk modulus K 0 parameters adapted to fit the data.

Noticing that the bulk modulus of a wide range of materials depends linearly on pressure, the EOS depends on the volume variation, the bulk modulus and its derivative. Usually, an EOS is associated to the volume variation under hydrostatic pressure. However, the consideration of the non hydrostatic contribution is poorly discussed in the literature but has a strong influence on the response of anisotropic materials. distance much higher than the in-plane intermolecular distances. Therefore, the initial black circular shape transforms into an elliptical shape, ensuring that the surrounding pressure field is purely hydrostatic. The Cauchy stress tensor applied on the single crystal then reads:

σ = ∆P • I, (6.2)
with ∆P the increment in pressure and I the second order identity tensor. Since in this case the sample behavior is controlled by pressure, the general strain-stress relationship is given by:

ε = ∆P S : I, (6.3)
where S is the fourth-order compliance elastic tensor, equal to the inverse of the elastic tensor (i.e. S = C -1 ), ε the second-order strain tensor and

(:) the inner product.

The bulk modulus could be deduced from S within the approximation of small strains (i.e. ∆P ≈ 0 and strain ε and deformation gradient tensor F related by F ≈ I + ε). The volume variation V V 0 = detF with V and V 0 the current and initial volumes, respectively, is given by:

V V 0 ≃ 1 + ε 1 + ε 2 + ε 3 , =⇒ ∆V V ≃ ε 1 + ε 2 + ε 3 , (6.4)
where only the first-order terms have been retained. In this configuration, the bulk modulus is given by:

K HP = V ∂P ∂V ≃ V ∆P ∆V . (6.5)
The combination of Equations 6.3 and 6.5 leads to the following expression of the bulk modulus:

K HP = 1 S 11 + S 22 + S 33 + 2S 12 + 2S 13 + 2S 23 . (6.6)
which depends on the sum of the upper left components of the compliance.

In the following, we choose to denote the EOS obtained through hydrostatic pressure "HP EOS".

Spherical deformation EOS

Considering the same reference state in Figure 6.1, the second type of loading investigated in this work is when a spherical deformation is applied to the sample. In this case, the deformation is fully controlled and the pressure is the observed variable. In the present case, the interlayer distance as well as the distance within the molecular layers undergo the same evolution, constrained by the application of the spherical deformation. Due to the very anisotropic properties of TATB single crystal, the resulting pressure field is not anymore hydrostatic. Therefore, the initial black circular shape in Figure 6.1 transforms into a smaller circular shape (i.e. homothetic transformation) with a non-hydrostatic surrounding pressure field.

Following the same line of reasoning as for the HP EOS presented before, we now assume that a spherical deformation is performed on a defect-free single crystal. Thus, the strain tensor reads:

ε = εI, (6.7)
and the Cauchy stress tensor is given by the following stress-strain relationship:

σ = εC : I. (6.8)
The combination of these expressions finally leads to the following bulk modulus:

K SD = C 11 + C 22 + C 33 + 2C 12 + 2C 13 + 2C 23 9 . (6.9)
where now the stiffness is summed up.

In the following, we choose to denote the EOS obtained through spherical deformation "SD EOS".

The two different loadings described above lead to different relations between pressure and volume for anisotropic materials. However, for isotropic materials, K SD = K HP and only one equation of state needs to be defined.

EOS experimental measurements

In order to study the equation of state of solids with diamond anvil cell (DAC) experiments, hydrostatic pressure is often expected. However, it may be difficult to achieve, due to both the type of sample and the pressure transmitting medium (PTM) used. This could lead to confusing results in the presence of anisotropic materials such as TATB. Ideally, if a single crystal is surrounded by a liquid PTM in a diamond anvil cell (DAC), hydrostatic pressure is guaranteed. However, if the PTM solidifies beyond a certain pressure, this leads to non-hydrostatic stress and non spherical deformation. Another method consists in loading the DAC with material powder in order to recolt X-ray patterns. If the powder is dispersed within a liquid PTM, hydrostatic conditions are fulfilled. However, if the PTM freezes or if grains are in contact in the DAC, non-hydrostatic stress is involved. The remaining of this chapter will present some preliminary results on TATB obtained with both molecular dynamics (MD) simulations and element free Galerkin method (EFG) with a comparison with experiments. The question to which this work aims to answer is the following:

When it comes to the mesoscopic modeling of polycrystalline materials for which the single crystal is very anisotropic, which equation of state is relevant ?

Molecular dynamics EOS computations

In the following, we describe two different techniques that have been used in our MD simulations of TATB single crystal, corresponding to the two equations of state presented above:

• The HP EOS: a hydrostatic pressure is applied

• The SD EOS: a spherical deformation is applied

Hydrostatic pressure EOS

We used the Parrinello-Rahman method in order to compute the HP EOS through MD simulations. This method (detailed in appendix A.

is equivalent to an equilibration in the ensemble NσT in which the Lagrangian equations are coupled with both a thermostat and a barostat, ensuring that σ ∝ I. This equilibration is performed through an equilibration of the cell parameters a, b, c and angles α β and γ, which implies:

σ = P cons I and F = V V 0 1/3 I, (6.10)
as discussed above.

Spherical deformation EOS

In order to compute the SD EOS however, we used the framework introduced in chapter 3 that allows the simulation box to follow any deformation path. A spherical deformation is applied to the MD simulation cell through a modification of the boundary conditions, and the overall stress tensor is computed during the trajectory. Thus, conversely to the computation of the HP EOS, the material undergoes the same amount of deformation in the three directions of the space, resulting in non-hydrostatic stress tensor due to the anisotropy of TATB single crystal. In a nutshell:

F = V V 0 1/3 I and σ = P I, (6.11)
as discussed above.

Results

The HP EOS was computed for pressures up to 50 GPa with a system size of 20a × 20b × 28c and is found to be in good agreement with the experimental data available in the literature [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF][START_REF] Olinger | The hydrostatic compression of explosives and detonation products to 10 gpa (100 kbars) and their calculated shock compression : Results for petn, tatb, co2 and h2o. 6[END_REF][START_REF] Plisson | Equation of state of 1,3,5-triamino-2,4,6trinitrobenzene up to 66 gpa[END_REF]. The simulation results are reported in Figure 6.2 with open grey triangles. Additionally, the Reuss modulus K HP has been computed from the compliance tensor presented in chapter 3 in order to fit the MD data through a Vinet EOS of which the only fitting parameter is K ′ R . The fitted curve (along with the values of K HP and K ′ R ) is plotted in Figure 6.2 with the red curve. We see that for pressures under 10 GPa the HP EOS fit deviates from the MD data. This is certainly due to the planar geometry of molecules considered in this work that produces a slightly higher compressibility of the interplanar distance, resulting in a C 33 elastic constant higher than other results in the literature [START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF][START_REF] Mathew | Generalized stacking fault energies in the basal plane of triclinic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene[END_REF]. However, we keep our value of K HP to be consistent with the complete elastic tensor that will be used in the following.

The SD EOS was computed in MD simulations with the same system size and for the same range of pressures. During the loading, a spherical deformation is prescribed to the TATB single crystal through MD simulations, starting from the same equilibrated crystal structure at (300 K, 0 GPa) as the HP EOS. The MD simulation results, reported in Figure 6.2 with open grey circles, can no longer be represented by a single EOS. Thus, we have chosen to split the MD data in three parts which are separately discussed below.

Interval V V 0 ∈ [1.0, 0 .88]
: SD EOS "compatible" evolution This interval corresponds to the beginning of the spherical deformation and no structural changes are observed. In this range of compression ratio, the deformation does not trigger any instability, despite an evolution of the stress components larger within the layers than perpendicular to it.

However, the evolution of the overall stress being quite "EOS compatible", a fit with a Vinet EOS is performed on this interval which can interpreted as a SD EOS. For this purpose, the modulus K SD has been computed from the elastic tensor of chapter 3 (Table 3.4) and the fit is performed by look-ing for an appropriate value of K ′ V . This SD EOS is represented with an orange line in Figure 6.2 as well as the values for K SD and K ′ V .

Interval V V 0 ∈ [0.88, 0.81]: Twinning-Buckling instability nucleation When V V 0 reaches a threshold value (i.e. ≈ 0.88), an instability is observed, leading to a volume decrease. This instability is similar to what was observed in chapter 5 with uniaxial compressions. Once the Twinning-Buckling (denoted TB in the following) is nucleated everywhere in the material (i.e. at V /V 0 = 0.81), the layered properties of the material are no longer present due to the strong misorientation provoked by the instability. Instead, a succession of chevrons with different orientations is observed and the elastic properties are expected to be more isotropic (e.g. the easy gliding in the basal plane before the instability can no longer be activated). This interval of "isotropization" however does not allow a fit to be performed due to its instable evolution. The resulting microstructure of such a spherical deformation is studied in the following.

Interval V

V 0 < 0.81: Layered to TB pseudo-phase transition At the beginning of this interval, the material has undergone a full TB instability. Because of both volume drop and an expected discontinuous evolution of the elastic properties (especially the shear components due to the large material rotations), this can be interpreted as a pseudo firstorder phase transition. The term "pseudo" is used since equilibration in the N σT ensemble should be investigated in order to confirm the volume discontinuity, signature of a first order phase transition. In the following, we prefer to term it "pseudo-phase transition" because this transition occurs in spherical deformation only and does not imply a change of the lattice but a new microstructure. Indeed, the local structure is still organized in molecular layers but the paving of the basal plane with chevrons of different orientations give different macroscopic properties to the TATB single crystal. This pseudo phase is named TB in the following.

Since for the TB phase the evolution of the overall stress is again "EOS compatible", we have decided to perform another fit while authorizing the optimization of both bulk modulus K T B and its pressure derivative K ′ T B . However, since no information is available on the volume at zero pressure, we have decided to keep it equal to the previous EOS. The corresponding curve is represented with the grey line in Figure 6.2 with its optimized parameters K T B and K ′ T B . The bulk modulus K T B is higher than both K HP and K SD but K ′ T B which relates to the high pressure material behavior is very close to K HP .

TB microstructure under spherical deformation

An additional MD simulation was performed in order to characterize the Twinning-Buckling pseudo-phase with a system size of 100a × 100b × 10c, equals to the size used in chapter 5. The system was strained instantaneously following a spherical deformation at V V 0 = 0.85. Then, an equilibration in the NVT ensemble was performed for 800 ps until the microstructure was equilibrated.

In chapter 5, we concluded that TATB single crystal under uniaxial compression systematically results in a Twinning-Buckling mechanism.

However, we did not study its behavior under biaxial or triaxial compression. We have found that even if a compression component is involved perpendicularly to the molecular layers, the Twinning-Buckling mechanism occurs, resulting in an equilibrated complex microstructure. The compressio ratio V V 0 = 0.85 is equivalent to ≈ 5% longitudinal strain in each direction, higher than the threshold around 3% estimated for uniaxial strain. The resulting microstructure after equilibration involves mainly two non-trivial deformation components of the Cauchy-Green strain tensor, that is C 13 and C 23 . In Figure 6 In order to characterize the microstructure and to link it with the twinning landscape calculated in chapter 5, we perform the analysis described in appendix B that consists in the processing of peaks of a density map d(C 13 , C 23 ) where the number of molecules with deformation of

(C 13 ± ∆ 2 , C 23 ± ∆ 2
) is given by d(C 13 , C 23 ) × ∆ 2 . The peaks locations as well as the associated proportion of molecules are reported in Figure 6.4. The total number of molecules associated with various deformation peaks is greater than 95%. Two main observations emerge:

• The two yellow circular dots are aligned with the strain-free minimum (C 13 ≈ C 23 ≈ 0) and centered around it with a deformation mainly in C 13 . This means that when interfaces have a [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] normal, pure buckling is involved. Indeed, in this direction, the elastic energy well is quadratic (see Figure 5.14 in the previous chapter), leading to the same amount of deformation but with opposite signs at each side of the interface (i.e. +∆C 13 and -∆C 13 ). This pure buckling mechanism concerns around 40% of the entire simulated material.

• The two white circular shapes however are only in C 23 deformation.

The larger shape, in the vicinity of the strain-free minimum (C 13 ≈ ize the microstructure fraction with positive and negative deformation during buckling and twinning. We note C T 23 the position of the twinning minimum. Then, ∆C 13 is the distance between the pure buckling deformation signatures and the strain-free energy minimum. Finally, ∆C 23 is the distance between the twinning deformation signature (small white shape) and the twinning minimum equal to the distance between the large white shape and the strain-free minimum. Since we have a pure compressive macroscopic deformation, the average shear deformation components are zero. For a compression direction along x we thus have:

< C

13 >= α∆C 13 + (1 -α)(-∆C 13 ) = 0, (6.12)

with α the material proportion with positive deformation. This expression leads to the trivial solution α = 0.5. This is consistent with the MD simulation results where the two yellow shapes in Figure 6.4 represents approximately 20% of the system each. Concerning a compression along the direction y, the two minima have to be taken in account. Thus we have:

< C 23 >= α(-∆C 23 ) + (1 -α)(C T 23 + ∆C 23 ) = 0, (6.13)
with α the material proportion with negative deformation, corresponding to the larger white shape in Figure 6.4. This equation leads to the final general solution:

α = C T 23 + ∆C 23 C T 23 + 2∆C 23 , (6.14)
which for C T 23 = 0 leads to the same result as before. In the case of TATB, C T 23 ≈ 0.5 and ∆C 23 ≈ 0.075. This leads to α ≈ 0.89. From the MD simulation the larger white shape corresponds to ≈ 90% of the overall white histograms, showing that the resulting microstructure proportion for buckling and twinning is consistent with this theoretical estimation.

We have shown that spherical deformation and hydrostatic compression of TATB single crystal lead to very different behavior. With hydrostatic pressure, TATB single crystal is very stable and no structural phase transition is observed. However, when a spherical deformation is applied, a very surprising Twinning-Buckling pseudo phase transition occurs and is expected to strongly modify the macroscopic elastic properties (at the microstructure level) of TATB due to the strong misorientations present in the microstructure. This may concern a wide range of materials with transerse isotropy. Additionally, this unusual mechanism could have happened during DAC experiments when polycrystals or powders were considered, leading to non hydrostatic pressure on the sample and thus potential misinterpretation.

TATB experimental results

To the best of our knowledge, only three experiments have been performed so far to determine the equation of state of TATB under hydrostatic compression.

Olinger and Cady [START_REF] Olinger | The hydrostatic compression of explosives and detonation products to 10 gpa (100 kbars) and their calculated shock compression : Results for petn, tatb, co2 and h2o. 6[END_REF] obtained X-ray powder diffraction patterns of TATB at different well-defined pressures up to ≈ 7 GPa. The explosive powder was placed in an annulus transparent to X-rays which was pressed with high uniaxial stress between Bridgman anvils. In their paper, they mention that the TATB powder was not packed but loosely inserted in the anvil cell, leading to hydrostatic compression conditions, in accordance with the HP EOS introduced above. During the determination of TATB cell parameters, the authors made two assumptions: the ratio a/b is kept constant with pressure and no relative displacements between two consecutive molecular layers is allowed. These assumptions lead to the necessity to measure only two independent parameters a and c. The data of Olinger and Cady is reported in Figure 6.5 as well as the two fits of HP EOS and SD EOS and will be discussed later.

Stevens et al. [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] performed powder X-ray diffraction experiments using a diamond anvil cell (DAC) and obtained the TATB unit cell volume for pressures up to 13 GPa. TATB crystals were ground with a mortar prior to loading. They used a Merrill-Bassett DAC with an Inconel gasket loaded with TATB powder, a small ruby chip and a 4 : 1 methanol/ethanol solution to ensure quasi-hydrostatic pressure conditions. However, they observe a freezing of the ethanol/methanol solution above 10 GPa, thus leading to non-hydrostatic conditions. Their results are reported in Figure 6.5.

Plisson et al. [START_REF] Plisson | Equation of state of 1,3,5-triamino-2,4,6trinitrobenzene up to 66 gpa[END_REF] determined the equation of state of TATB up to 66 GPa using synchrotron X-ray diffraction experiments in two different DAC. The first DAC cell was loaded with a TATB powder without pressure transmitting medium and intended to obtain diffraction patterns at ambient pressure to compare with literature. The second DAC cell was loaded with a small TATB sample that consisted in a few single crystals surrounded by neon as pressure transmitting medium. The obtained values are reported in Figure 6.5 and cover a wider range of pressure than the two previous experiments.

It is clear that among the three expirements, the one performed by

Mesoscopic modeling of a TATB polycrystal by non linear elasticity

In the following, we propose a first modeling of the TATB single and polycrystal mechanical response by using non linear elasticity in a Lagrangian formulation. The main goal of this section concerns the modeling of an elastic instability under spherical deformation of TATB single and polycrystal apparented to a pseudo-phase transition, observed through molecular dynamics simulations.

During this pseudo-phase transition, the twinning energy landscape has an influence only on the resulting microstructure and not on the nucleation process, due to pure elastic considerations. Since we focus on the nucleation of such an instability and not on the resulting microstructure, we decided not to take into account the twinning energy landscape as a first approximation. Thus, we use the term Elastic Buckling simulations (EB simulations) to denote the results obtained within the Lagrangian continuum mechanics framework.

We need first to derive the stiffness in the context of large strain.

Then, the model will be validated on the behavior of the TATB single crystal by direct comparison between MD and EB simulation. Finally, the model will be applied to a polycrystalline aggregate of TATB.

EOS derived stiffness tensor

The aim of this section is to deduce the stiffness tensor, function of the volume variation V /V 0 , from a EOS [START_REF] Denoual | Personnal Communication[END_REF] in the context of large strain formulation. In order to do this, we need:

• The stiffness tensor C 0 = C 0 (V /V 0 ) = C 0 (1) at ambient conditions.

• The bulk modulus B 0 at ambient conditions obtained from C 0 .

• The spherical deformation EOS with the relation

P EOS = f (V /V 0 ).
with f a regular function usually tabulated from the MD or experimental data. In the present case, we consider a spherical deformation (i.e. F = By imposing F = (V /V 0 ) 1/3 I (spherical deformation), we obtain the following relation for P :

P = F • D : E + e 2 2 I : D ′ : I d(V /V 0 ) dF (6.20)
with e = 1 2 ((V /V 0 ) 2/3 -1) and E = eI. Additionally, we have:

(F • D : E) mn = e √ 2e + 1D mnkk d(V /V 0 ) dF mn = V /V 0 • F -T = (2e + 1)δ mn (6.21)
leading to the following expression for P :

P mn = e √ 2e + 1D mnkk + e 2 (2e + 1) 2 D ′ mnkk δ mn . (6.22)
Then, with the usual relation between the Cauchy stress tensor σ and

P σ = 1 V /V 0 F P T , (6.23)
and using in indicial notation, the latter expression reads:

σ pm = 1 2e + 1 δ pn P mn , = e √ 2e + 1 D mpkk + e 2 2 D ′ mpkk δ mp . (6.24)
By imposing the relation P EOS = 1 3 trace(σ) = 1 3 σ ii and using the expressions in Equation 6.19, we obtain:

P EOS = e 3 √ 2e + 1 D ijkk + e 2 6 D ′ ijkk δ ii , = 3e √ 2e + 1 B + 9e 2 2 B ′ , (6.25)
which is equivalent to the following first order differential equation:

B ′ (V /V 0 ) + 2 3e √ 2e + 1 B(V /V 0 ) = 2 9e 2 P EOS (V /V 0 ). (6.26)
Finally, the volume dependent stiffness tensor D reads:

D(V /V 0 ) = B(V /V 0 ) B 0 C 0 (6.27)

Validation on the TATB single crystal

In this section, we aim at reproducing the behavior of a TATB single crystal subjected to a spherical deformation which was characterized through MD simulations in section 6.3.3. We are interested in reproducing the spherical deformation equation of state as well as the buckling instability nucleation threshold. However, in the following, we do not have taken in account the twinning energy landscape. Thus, we expect to have some differences on the resulting microstructure.

An element-free Galerkin least-squares formulation [START_REF] Belytschko | Element-free galerkin methods[END_REF] The simulated material is a cube of TATB single crystal containing 64000 elements with 3D periodic boundary conditions of fictive total length L = 3.9 mm. The final time t f of the dynamic continuum mechanics simulation used in this work is related to the time t c needed by acoustic waves to go through the entire sample. In the present case, quasistatic loading conditions are ensured by assuming t f = 100t c , i.e. the acoustic waves run over 100 times through the sample during the en-tire simulation. Thus, t c = L/c L with c L the wave velocity in TATB, i.e., c L = K SD /ρ 0 = 3043 m s -1 , and t f ≈ 128 µs.

A two-steps loading is then proceeded on the initial TATB single crystal:

1. A spherical deformation up to V V 0 = 0.85 is achieved during the interval of time t 0 < t < 0.1t f . Remark: This is exactly equivalent to the method used to measure the spherical deformation equation of state through MD simulations up to V V 0 = 0.85. In Figure 6.6 are reported the stress/strain histograms (colored points with size dependent on a volume fraction fo the system) of the EB simulation described above as well as the average stress/strain values (grey empty squares). We can see that up to V V 0 ≈ 0.88 both histograms and average values follow the spherical deformation EOS (SD EOS), confirming that the computation of pressure-dependent elastic constants is well defined.

During this interval, the TATB single crystal is uniformely strained and no elastic instability is generated. However, when V V 0 < 0.88, the histograms become very dispersed and the average values go between the SD EOS and the HP EOS. In that case, as the specimen becomes "polycrystalline", SD and HP EOS can be considered as Voigt upper bound and Reuss lower bound, as they involve average of C and S, see equations 6.9 and 6.6. This EB simulation reproduces correctly the volume-pressure evolution computed through MD simulations. Furthermore, the elastic buckling instability is triggered at a compression ratio similar to that of MD which shows that the non-linear elasticity framework used is consistent with the SD EOS computed in MD.

At t = 0.1t f , the volume is kept constant with V V 0 = 0.85 and the microstructure relaxed up to t = t f . Its evolution is shown in Figure 6.6 under the label "EB Simulation" where the elements E 13 of the Green-Lagrange strain tensor has been chosen for the color coding. Under the label "MD Simulation", the microstructure after 800 ps of relaxation is shown, with the same color coding.

When reasoning in terms of acoustic waves, the resulting microstructure of the MD simulation presented in Figure 6.6 corresponds to a time after which acoustic waves have run approximately 25 times through the sample 2 . Since during the EB simulation, the final time equals 100 runs of the acoustic waves through the material, the equivalent microstructure that should be compared with MD simulation is the one at t = 0.25t f . For the EB simulation, the entire single crystal volume is affected by a nonzero deformation in E 13 . Conversely, ≈ 50% of the single crystal volume in the MD simulation is affected by zero deformation E 13 which concerns zones affected by a deformation E 23 corresponding to the twinning deformation. This highlights the differences between the activated mechanisms: during the EB simulation, a buckled microstructure is present everywhere in the single crystal volume whereas a combination of buckling and twinning is observed in MD. As announced before, the twinning energy landscape was not taken in account in this study. Thus, the fully buckled microstructure with no twinning obtained with the EB simulation is consistent with what was expected.

The two simulations exhibit similar results. The general behavior is the same with the two methods and the nucleation of the elastic buckling instability is provoked at the same compression ratio. Furthermore, the resulting microstructure obtained with both simulations can be compared since acoustic waves have made the same number of round trips through both samples. Since the twinning energy landscape was not considered here, we show that the EB simulation produces a fully buckled microstructure whereas a coupling between twinning and buckling was observed in MD simulations. Thus, the EB simulations are in qualitative agreement with molecular dynamics, allowing to investigate the mechanical behavior a small TATB polycrystal under macroscopic spherical deformation.

Simulations of a TATB polycrystal

In order to generate the 3D periodic polycrystal, we first generate a list of 10 randomly distributed crystal with different orientations. Then, 10 seeds are randomly distributed within the cube and correspond to the center of each different future grain. Finally, a Voronoi tessellation is performed to generate the polycrystal from the seeds with specific orientations. Quasi-static loading conditions are ensured as for the single crystal simulation and the same two-steps loading is adopted up to V V 0 = 0.80 (i.e. spherical deformation for t 0 < t < 0.1t f , equilibration for 0.1t f < t < t f with the same time and length values).

In Figure 6.7 are reported the stress/strain histograms for the polycrystalline TATB elastic buckling (EB) simulation along with their average values. The color coding used is the same one as in Figure 6 and we can see that for the step 5 (t = 0.09t f ), the histogram is very spread and lies between the SD and HP EOS.

Views of the microstructure corresponding to the different steps reported in the graphs are presented in Figure 6.7. We can see that the first signs of an elastic buckling instability within grains can be detected at t = 0.05t f which corresponds to V V 0 = 0.90. This value is similar to the threshold for the TATB single crystal spherical deformation simulations using both molecular dynamics and continuum mechanics techniques (< 0.88). Thus, the spreading of the pressure-volume histograms are due to the elastic interactions between grains, and, even if the macroscopic applied deformation is spherical, each grain experiences a non spherical deformation due to the random orientations of its neighbors. As the deformation continues, we can see that the buckling grows and that every grain is buckled. Furthermore, we think that the interfaces are oriented perpendicularly to the basal plane from what was observed in MD simulations. The color coding for the microstructure is based on the element E 13 of the Green-Lagrange strain tensor. Similarly to the single crystal, all the crystals are entirely buckled, which is consistent with the absence of a twinning energy landscape.

The simulation of the TATB polycrystal shows that the SD EOS computed on the single crystal is not reproduced by the polycrystal behavior for which the pressure-volume relation deviates early from the EOS curve.

However, the nucleation of the elastic buckling instability is observed at a similar compression ratio than for the single crystal which shows that non-hydrostatic stresses are generated in the entire polycrystal, due to the anisotropic elastic interactions between grains.

Conclusion

In this chapter, we have presented the behaviors for two loadings conditions, i.e. when a TATB single crystal is subjected to hydrostatic pressure conditions and to spherical deformation.

We performed molecular dynamics simulations in order to show these differences for the TATB single crystal. An equation of state under hydrostatic pressure up to 60 GPa was computed using Parrinello-Rahman simulations. The TATB single crystal is found to be very stable while maintaining its anisotropic layered structure. Additionally, thanks to the prescribed deformation paths methodology introduced previously, we performed spherical deformation on a TATB single crystal, starting from the same structure at classical ambient conditions. This deformation exhibits an "EOS compatible" evolution for 1.0 > V /V 0 > 0.88 while for 0.88 > V /V 0 > 0.81 Twinning-Buckling is observed. Once this instability stabilizes, an evolution compatible with an EOS is observed for V /V 0 < 0.81. The two different equations of state are denoted HP EOS for hydrostatic pressure and SD EOS for spherical deformation.

For a spherical deformation with V /V 0 = 0.85 the resulting microstructure is a coupling between pure twinning and pure buckling. This result is very surprising because it shows that when spherical deformation is applied to very anisotropic materials, a pseudo-phase transition can be observed, which is impossible to trigger under hydrostatic pressure.

The equation of state computed with molecular dynamics was then used to set up a non-linear elasticity modeling in a Lagrangian framework.

Then, a spherical deformation equals to the one applied in MD simulations was applied to a TATB single crystal. A good agreement was found between the continuum mechanics simulation and the MD data. The buckling instability is reproduced at the same threshold of V /V 0 ≈ 0.88.

However, in MD simulations a coupling between twinning and buckling is present, that cannot be reproduced as the twinning energy landscape was not taken in account in the Lagrangian framework which focused only on the generation of the buckling instability.

The qualitative agreement between non linear elasticity and MD allowed us to consider the polycrystal case using the same approach. A polycrystal containing 10 randomly oriented grains of TATB single crystal was created and a macroscopic spherical deformation applied. Due to the presence of multiple grains with different orientations, non-hydrostatic pressure is generated and the elastic buckling instability is observed at the same threshold of V /V 0 ≈ 0.90. However, the stress/strain local states are found to be much more dispersed but remain between the two bounds defined by the hydrostatic pressure and spherical deformation equation of states. Finally, the resulting buckled microstructure is present within every grain showing that none of them was subjected to pure hydrostatic pressure.

These different results represent an important step toward the modeling of TATB single and polycrystal mechanical behavior. A few comments have nonetheless to be made. Taking the spherical deformation EOS computed from MD simulations as an input in the continuum mechanics framework lead to a very good reproduction of the pressurevolume relation and a similar microstructure. The differences between SD and HP EOS should vanish for a polycrystal with randomly oriented crystals as the effective behavior to consider in the EOS becomes isotropic.

Adapted homogeneization methods are then needed to compute the local stress-strain behavior inside each grain, and its fluctuation, in order to predict the final buckled and twinned microstructure. However, a pseudophase transition such as buckling of each grain in the polycrystal shows the complexity of its behavior and demonstrates the necessity (and the challenge) to consider homogeneization methods compatible with this "pseudo" phase change in order to compute a representative equation of state. Dans ce chapitre, nous nous focalisons sur les matériaux très anisotropes, en prenant le TATB comme exemple, et sur leur comportement sous contrainte hydrostatique et non-hydrostatique, menant à des différences remarquables.

Résumé substantiel en franc ¸ais

Équations d'état en pression hydrostatique et en déformation sphérique pour les matériaux très anisotropes

La caractérisation de la densification des matériaux sous pression a fait l'objet de nombreuses études, donnant lieu à des fonctions empiriques telle que l'équation d'état de Vinet [START_REF] Vinet | A universal equation of state for solids[END_REF] :

P = 3K 0 x 2 (1 -x) exp 3 2 (K ′ 0 -1)(1 -x) (6.28) avec x = V V 0 1/3
σ = P cons I et F = V V 0 1/3 I.
En utilisant la méthode présentée dans le chapitre 3 permettant à la boîte de simulation de suivre un chemin de déformation quelconque, on peut appliquer une déformation sphérique au matériau simulé menant à une contrainte non-hydrostatique. On a donc dans ce cas là : Classical molecular dynamics technique was used in order to simulate the TATB single crystal at various thermodynamic conditions (temperature and pressure) and under different loadings, in order to get some insight into physical and mechanical properties. Specific diagnostics, in particular mechanical tools, were implemented in MD simulations so as to inform continuum mechanics simulations (i.e. at the mesoscopic scale), which represents an inevitable step if one wants to build a representative constitutive law for polycrystalline TATB. These tools deepen our understanding of the TATB behavior at the microscopic scale and happen to be essential for the comparison between microscopic (nm) and mesoscopic (µm) scales, that share the same diagnosis.

F = V V 0 1/3 I et σ = P I.

Les résultats présentés en

D(V /V 0 ) = B(V /V 0 ) B 0 C 0 ( 6
An all-atom potential from the literature was used to simulate the equilibrium states of TATB single crystal through molecular dynamics simulations using rigid molecules in chapter 3. The cell parameters were computed using the Parrinello-Rahman method up to 50 GPa and were found to be in good agreement with the experiments performed at CEA-DAM.

Starting from the equilibrium structure at any (T, P ) conditions, a method to compute the elastic constants in temperature was introduced and the results were found to be in satisfactory agreement with the literature. A similar method was also presented and allowed for the computation of TATB basal planes P 1 and P 2 γ-surfaces in temperature and pressure.

These preliminary results validate the rigid molecules approximation for computing mechanical properties.

Within layered materials, dislocations are known to easily nucleate in the basal plane (i.e. between the atomic/molecular planes). Thus, the dislocation core structures for [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] and [010] Burgers vectors in P 1 and P 2 planes were investigated using the generalized Peierls-Nabarro-Galerkin (PNG) model in chapter 4, which represented a first attempt to bridge scales. Due to the strong elastic anisotropy of the TATB, the calculation of core structures of such dislocations is out of reach in MD.

Indeed, the characterization of such dislocations would require large system size (on the order of 10 7 unit cells), thus making a direct approach by MD simulation hardly possible. Instead, we have used MD to compute the necessary ingredients for a mesoscopic simulation using the PNG method. Elastic constants as well as γ-surfaces at 300 K and for pressure up to 10 GPa were used as inputs in the PNG model in order to compute the corresponding dislocation core structures. Dislocations (edge and screw) in P 1 and P 2 planes were found to split into up to four partials and exhibit a very complex core, the distance between the two most distant partials being around a hundred Burgers vectors. Furthermore, pressure was shown to have no significant effect on the dislocation core structures, due to the high stability of TATB single crystal under hydrostatic pressure.

Thus, dislocations in TATB basal planes can not be considered as discrete lines but rather as large stacking fault ribbons for the range of pressure studied. The ability of TATB to glide easily on its basal plane raises the question of the computation of a viscosity parameter for gliding in these directions in order to model this behavior at the mesoscopic scale.

In order to impose deformations spanning all relevant directions to single out the mechanisms of deformation, we have implemented a method that allows to prescribe any deformation to the MD simulation supercell following prescribed paths, i.e. directional pure shear, uniaxial and/or triaxial compressions, etc., in a similar way to what is used in continuum mechanics Lagrangian codes. Furthermore, since large deformations of materials may lead to the nucleation of defects such as dislocations and twinning that involve strong local perturbations and rigid body rotations, we wanted to be able to identify and characterize these processes. To this aim, we have implemented in the MD code a method that allows to compute at each particle location the local deformation gradient tensor with respect to a reference state. From this deformation gradient tensor, a polar decomposition can be done in order to access to the entire strain history of the material, i.e. planes that have undergone dislocation gliding or zones with strong rotations, a typical signature of twins.

The prescribed deformation paths method was used to apply directional pure shear to a TATB single crystal in chapter 5, using a compression coupled with an equivalent transverse expansion. A yield von Mises stress surface was computed and led to the identification of three distincts behaviors. When pure shear is applied to the basal plane, easy gliding is activated and no other deformation mechanism than stacking fault ribbons is observed.

However, when pure shear is applied in the direction of non-basal planes, MD simulations predict nucleation of dislocations. Four independent slip systems were identified. These dislocations are nucleated through a hydrogen-bond breaking process and have complex cores with a local dilatancy up to 10 -15%. Dislocations are known to be the main mechanism that controls plastic activity in metallic alloys and minerals.

However, in the aim of building a constitutive law that would involve a crystal plasticity model, many difficulties have to be overcome, e.g., obtaining data on dislocations velocities, interaction modes, core structures and activation stresses etc. Furthermore, the majority of such properties could not be studied with MD simulations, due to the need of very large systems, to be able to simulate long range elastic interactions between dislocations for example. The dependence of dislocation nucleation with pressure was studied, resulting in larger nucleation stresses for increasing hydrostatic pressure. Then, the question of considering non-basal dislocations within consitutive laws for the behavior under shock loading remains open, regarding the stress needed for their activation.

When the compression direction lays in the TATB basal plane, a surprising mechanism is involved: when the longitudinal strain reaches ≈ 3%, an elastic bukling instability is observed. Such a mechanism, which consists in a succession of chevrons with positive and negative deformation was found to be very dependent on the loading direction. The interfaces between chevrons are always normal to the loading direction and the lower the strain rate is, the larger the chevrons are. Such a deformation can be locally seen as a composition of pure shear followed by a rotation.

In order to understand the resulting microstructure, we have computed a homogeneous shear energy surface by applying a homogeneous shear within the basal plane of TATB single crystal. This computation results in the identification of a twin deformation path of the TATB lattice. In order to find a correlation between this energy landscape and the microstructure observed when buckling is triggered, we extracted the deformation distribution from the MD simulation by using the local deformation gradient tensor. The two non-trivial components of the deformation involved during buckling, reported on the energy landscape, show that the buckling behavior strongly depends on this landscape. Indeed, the deformation states are found to be localized in the vicinity of the two minima present on the landscape. Thus, when compression is applied along the [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] direction, pure buckling is involved. Conversely, the more the compression direction is oriented perpendicularly to the [START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] direction, the more twinning is activated.

This new mechanism was denoted as Twinning-Buckling mechanism. Thanks to previous work on the buckling mechanism of layered materials and block copolymers, we were able to study the nucleation property of buckling in TATB in linear elasticity. Considering buckling as a succession of shear and rotation, a link between the local deformation components and the macroscopic ones is drawn, corresponding to pure compression with a small transverse perturbation. We found a theoretical threshold around 4% for this instability, of the same order than the one observed during MD simulations. The fact that the nucleation of this buckling instability comes from purely elastic consideration is very interesting and surprising at the same time: whenever a compression with a longitudinal strain above 3% is applied within the TATB molecular layers, the Twinning-Buckling mechanism is activated and strongly modifies the local microstructure. This mechanism could prevail over all the other mechanisms presented early in the first steps of the mechanical response of TATB.

As a perspective, a constitutive law that would take in account all the identified deformation mechanisms would be very interesting to construct. However, since we do not have enough informations on the nonbasal dislocations at the moment, classical crystal plasticity can not be envisaged. A solution would be to consider a plastic flow law such as Norton-Bingham in order to model the plastic activity in TATB single crystal as well as its pressure dependence (accounting for the local dilatancy involved in the nucleation process). The nucleation of the buckling instability being purely elastic, it is well reproduced with non linear elasticity as shown in the last chapter of this thesis and discussed below. Finally, we showed that the microstructure after buckling strongly depends on the Twinning energy landscape. A solution to make possible the twinning phase transition during buckling would be to use phase field reaction pathways as done by Denoual et al. [START_REF] Denoual | Phase-field reaction-pathway kinetics of martensitic transformations in a model fe 3 Ni alloy[END_REF] for martensitic transformation in a model Fe 3 Ni alloy and by Denoual [START_REF] Denoual | A phase field approach with a reaction pathways-based potential to model reconstructive martensitic transformations with a large number of variants[END_REF] and Vattré [START_REF] Vattré | Polymorphism of iron at high pressure: A 3d phase-field model for displacive transitions with finite elastoplastic deformations[END_REF] for the iron behavior under pressure.

In chapter 6, we focused on the computation of an equation of state for very anisotropic polycrystal mechanical behavior with an application to TATB. Two equations of state were considered, leading to very different results. When hydrostatic pressure is applied on the single crystal, TATB cell parameters undergo different evolutions, resulting in a non spherical deformation. Conversely, when a spherical deformation is applied, each cell parameter undergoes the same amount of deformation and nonhydrostatic stress is involved. The hydrostatic pressure equation of state (HP EOS) was computed using Parrinello-Rahman MD simulations while the spherical deformation (SD EOS) was computed using the prescribed deformation paths context starting from the same equilibrated TATB single crystal. At a compression ratio of ≈ 0.9, a pseudo-phase transition corresponding to the Twinning-Buckling mechanism was triggered during the SD EOS computation but was not observed during the HP EOS computation, due to the high stability of TATB under hydrostatic pressure. This implies that when very anisotropic materials behavior under pressure is investigated, a pseudo-phase transition can be observed under spherical deformation but is completely absent with hydrostatic pressure. This shows the importance of non hydrostatic stress for phase transition in very anisotropic materials.

The behavior of TATB single crystal under spherical deformation was then simulated using large strain non linear elasticity along with the SD EOS computed with MD simulations as an input. A spherical deformation was prescribed up to V V 0 = 0.85, and imposed constant after. The volume-pressure relation was in qualitative agreement with the SD EOS obtained with MD simulation and the buckling instability reproduced at the same compression ratio. After the nucleation of the buckling instability, the average volume-stress deviates from the SD EOS and stay between the HP EOS and SD EOS. However, since twinning was not taken into account, the resulting microstructure was only composed of buckled parts, in opposition to the microstructure obtained with MD containing twinning as well. The comparison was made possible by the computation of the local deformation gradient tensor, in MD simulations, and similar amplitudes of deformation were observed.

Due to the qualitative agreement between MD simulations and non linear elasticity buckling (EB) simulations, we chose to simulate a small TATB polycrystal within the same framework. Ten grains with random orientations were generated using a Voronoï tessellation to create a polycrystal. Then, a macroscopic spherical deformation was applied. Due to the random orientation of each grain, elastic interactions generate nonhydrostatic stresses and the buckling instability is triggered at the same compression ratio than simulations on the single crystal. However, the average volume-stress distribution deviates very early from the SD EOS and stay between the HP and SD EOS as well. This simulation led to the presence of a fully buckled microstructure in each grain, showing that the macroscopic spherical deformation implied non-hydrostatic stresses in the entire polycrystal. These different results show the complexity of the polycrystal behavior for TATB. The necessity to consider homogeneization methods that would take in account this pseudo phase transition represents a real challenge and should lead to very interesting results in the future.

Conclusion et perspectives en franc ¸ais

Au cours de cette thèse, une approche multiéchelle a été utilisée q 2 1 + q 2 4 -1 2 q 1 q 2 + q 3 q 4 q 1 q 3 -q 2 q 4 q 1 q 2 -q 3 q 4 q 2 2 + q 2 4 -1 2 q 2 q 3 + q 1 q 4 q 1 q 3 + q 2 q 4 q 2 q 3 -q 1 q 4 q 2 3 + q 2 4 -

1 2     (A.2)
With this definition, the quaternions of each molecule verifies easily the equations of motion which no longer include singularities.

A.2 Parrinello-Rahman Simulations

Parrinello-Rahman simulations [START_REF] Parrinello | Crystal structure and pair potentials: A molecular-dynamics study[END_REF][START_REF] Parrinello | Polymorphic transitions in single crystals: A new molecular dynamics method[END_REF] are very usefull to find the equilibrium structure of crystals at different conditions of temperature and pressure, allowing for example to find the different polymorphs of certain crystals. We briefly describe the method here since it has been used during this work to find TATB cell at various temperature and pressure conditions.

A crystal primitive cell can have arbitrary shape and volume completely described by its three lattice vectors {a, b, c}, of different lengths and not necessarly orthogonal with each other in the most general, triclinic case. A 3 × 3 matrix h is constructed (sometimes denoted "lattice matrix"), its columns being equal to the components of a, b and c. Cell volume Ω, position of particle r i and square of distance between i and j particles r 2 ij can then be defined as: Unfortunately, due to the perturbed trajectories of MD simulations and the large deformation spectrum that is spanned by the overall system, the resulting density map is most of the time very noisy. This is the case of our example for which the density map of the initial data is shown in In order to overcome this problem, a smoothing is performed on the initial data, through a FFT convolution with a gaussian kernel. This results in a very smooth density map (see in Figure B.1) on which various operations can be easily done such as identify the exact positions of different significative density maxima (four in number). This procedure is described below.

Ω = ||h|| = a • (b ∧ c), r i = hs i = ξ i a + η i b + χ i c,

B.2 Maxima Locations Identification

Since we now have a density map which allows us to identify various distinct maxima, the aim is to compute their exact position, i.e. a couple C 13 , C 23 . The first step consists to skim the density landscape from the highest values to the lowest ones. This is equivalent to bringing down a "sheet" from altitude z = max(d(C 13 , C 23 )) to z = 0. Doing this, the "sheet" will interact with the various peaks, starting with the one of maximum density. A new zone is identified each time an interaction occurs (see

B.3 Association Between Positions and Density

A simple maximum extraction is then performed for each zone and the positions of the different maxima can be reported onto the initial density map, represented as star shapes in

Multiscale Modeling of Energetic Materials Behavior ABSTRACT:

The construction of mesoscopic (micrometer scale) constitutive laws in materials science is studied for a long time. However, the constant progress in high performance computing changes the perspectives. Indeed, constitutive laws now aim at explicitly take into account the microstructure and its underlying physics at the atomic scale, for which simulation techniques prove to be very accurate but definitely expensive. The multiscale approach is therefore perfectly adapted to such a challenge and the dialogue between scales necessary. In this thesis, the mechanical behavior of the energetic material TATB in temperature and pressure is investigated using molecular dynamics simulations in order to understand the microscopic deformation mechanisms responsible for plastic activity. The local computation of mechanical variables was developed in atomistic simulations, allowing the dialogue with continuum mechanical methods. Additionally, prescribed deformation paths were coupled with molecular dynamics, allowing to reveal the plasticity mechanism of TATB single crystal. Nucleation of complex dislocation structures with intrinsic dilatancy, twinning transition pathway and a twinning-buckling pseudo phase transition are three distinct behaviors triggered for different loading directions. Then, mesoscopic simulations inferred by atomic scale observations aim at reproducing the twinning-buckling pseudo-phase transition under tri-axial compression using a Lagrangian code. The comparison between both simulation techniques is made possible thanks to the mechanical tools that have been implemented in the molecular dynamics code. Finally, polycrystalline TATB is simulated with non linear elasticity and we demonstrate the necessity to consider an equation of state compatible with this pseudo phase transition, which has a strong influence on the polycristal behavior.

  Cette thèse a été réalisée dans le cadre d'un partenariat entre le CEA-DAM pour l'encadrement et les Arts et Métiers ParisTech pour la direction. L'encadrement a été effectué par M. Christophe DENOUAL et M. Jean-Bernard MAILLET, tous deux Ingénieurs de recherche au CEA-DAM. Côté Arts et Métiers, M. Olivier CASTELNAU, Directeur de recherche CNRS et Mme. Katell DERRIEN, Maître de conférences du laboratoire Procédés et Ingénierie en Mécanique et Matériaux ont endossé les rôles de directeur et co-directrice de thèse. Sans l'encadrement dont j'ai bénéficié, le travail présenté dans ce manuscrit n'aurait pas été possible et je leur suis particulièrement reconnaissant.

  rais remercier M. Nicolas PINEAU et Mme. Claire LEMARCHAND pour leur aide lorsqu'il a fallu mettre les mains dans le code de dynamique moléculaire mais aussi pour leur bonne humeur constante, qui a été un vrai plus ! Je remercie également M. Ronan MADEC et M. Yves-Patrick PELLEGRINI pour nos multiples interactions autour de ma thématique qui ont donné (et donneront certainement) lieu à de belles discussions. Ne pouvant pas citer tout le monde, je tiens à remercier l'ensemble du personnel du service pour tous les échanges qui ont pu avoir lieu lors des séminaires ou encore des pauses café. Une mention spéciale à Mme. Sandra BOULLIER dont l'aide et la compréhension ont fortement facilité l'aspect administratif. Avant de passer à des remerciements moins formels, j'aimerais remercier M. Franc ¸ois JOLLET et Mme. Catherine CHERFILS pour leur suivi de ces dernières années et pour la confiance qu'il m'ont faite dans l'optique de poursuivre l'aventure au sein de leur équipe. Place à l'ensemble des stagiaires/doctorants/post-docs et autres personnes avec qui j'ai passé d'excellents moments au cours de ces trois dernières années, lors des multiples repas, pauses cafés et j'en passe. Je les remercie particulièrement de m'avoir attendu chaque midi pour mes retards répétés, notamment David, Raphaël et Giovanni. Cette thèse n'aurait pas été aussi agréable sans leur présence et pour c ¸a, un gros merci à : Ahmed, Emmanuel, Gérôme, Jean-Charles, Jordan, David, Lucas, Luc, Aude, Richard, Xavier, Sami, Grégoire, Jean-Baptiste et Théo. Une mention spéciale à Estelle, Jean-Baptiste et Thibaud, camarades de thèse à qui je souhaite plein de belles choses pour la suite, aux deux terreurs Raphaël et Giovanni à qui je souhaite une belle dernière année et à Luis à qui je souhaite une belle thèse ! vi Je tiens bien évidemment à remercier mes amis "en dehors" du CEA qui ont toujours été là pour moi et qui m'ont aidé à me changer les idées lors des périodes difficiles. Merci aux copains du Sud-Ouest pour ces beaux moments passés ensemble lorsque j'avais l'occasion de redescendre, notamment au Heugais et à la meute ! Merci également aux camarades de l'ENSTA et amis, pour tous ces beaux moments passés ensemble lors de vacances et week-ends, je pense notamment à la team Kikou et à la famille ENSTA (presque) du Sud-Ouest. Bien évidemment, merci aux copains de l'équipe Radio Déchets qui ont rendu ces trois dernières années tellement agréables grâce à des repas plus ou moins gastronomiques, de belles soirées et quelques voyages culturels bien sympathiques ! Enfin, je remercie les membres la Team BB, partenaires dans la douleur mais surtout amis pour toutes ces heures de "Light weight" passés au Crossfit et ailleurs ! Et puis je remercie toutes les personnes que j'ai pu oublier, collègues thésards de Saclay, copains d'école, copains d'enfance. Et pour terminer, je voudrais remercier du fond du coeur mes parents pour toutes les choses qu'ils m'ont enseignées toutes ces années et sans qui je n'aurais jamais été capable d'atteindre mes objectifs. Je les remercie particulièrement de m'avoir toujours encouragé lors de mes choix d'orientation et de n'avoir jamais douté de moi. Leur soutien, leur accompagnement et leur aide ont été des éléments moteurs lors de ces dernières années. Je suis également très reconnaissant envers mes deux soeurs avec qui j'ai passé d'excellents moments en famille pour le courage qu'elles m'ont donné et pour leur soutien inconditionnel toutes ces années ! A votre tour maintenant

Figure 1 . 1 :

 11 Figure 1.1: TATB unit cell with lattice vectors a, b, c and cell angles α, β and γ. (a) 3D representation. (b) x+ view. (c) z+ view. (d) y+ view.

Figure 1 . 2 :

 12 Figure 1.2: Typical SEM images of different TATB samples with a ×500 (up) and a ×100 (down) magnification. Image after Guerain et al. [47]

  constructed following TraPPE potentials (Transferable Potentials for Phase Equilibria) of aniline and nitrobenzene molecules. On the opposite, Gee et al. [44] developed a potential for fully flexible molecules, allowing any kind of distortion. Finally, Bedrov et al. [11] introduced a FF (with both polarizable and non-polarizable versions) which consider bonds fixed at their equilibrium values. Since the non-polarizable version of the latter has been used in this work, it is detailed in section 3.1.2. The above enumerated FFs are all constructed in order to reproduce TATB structural properties such as cell parameters but some are parametrized only on quantum-chemistry calculations and other on additional experimental data. Nevertheless, they were (more or less) all fitted to reproduce quantum-chemistry calculations for the electrostatic potential of an isolated TATB molecule using partial atomic charges. The main differences in the construction of the above FFs were in the formulation used to model interatomic interactions. Gee et al. and Rai et al. both use Lennard-Jones functions for the repulsion and dispersion terms whereas Bedrov et al. employ Buckingham (exp -6) functions. The latter is described in detail in chapter 3. Due to the lack of experiments concerning thermo-mechanical properties of TATB single crystals, classical MD simulations have been widely

Figure 2 . 1 :

 21 Figure 2.1: TATB-based PBX microstructure. Optical microscopy micrograph obtained through reflected polarized non analyzed light. Image after Bruneton et al. unpublished.

Figure 3 . 1 :

 31 Figure 3.1: Pressure dependence of TATB cell lengths (a), cell angles (b) at 300 K and equation of state (c) from MD simulations compared with experiment and density functional theory data from Pineau et al.

Figure 3 . 2 :

 32 Figure 3.2: Evolution with pressure of the nine orthotropics stiffness coefficients of TATB single crystal at ambient temperature.

Figure 3 . 3 :

 33 Figure 3.3: Plots of elastic-energy vs. strain curves. Black dots: simulations data from multiple states. Continuous lines: final fitted curves. (a), (b), (c) and (d) correspond to C 11 , C 33 , C 44 and C 12 elastic constants, respectively.

Figure 3 .Figure 3 . 4 :

 334 Figure 3.4 and I is the identity tensor. Two types of loadings are used in the following: pure shear strain is obtained by taking β = -α/(1 + α) (implying that detF = 1), and uniaxial loading compression with β = 0.

  the dislocation line, but not the dislocation slip plane. For defects such as twins, the local perturbation of crystal order only gives information about the location of twins' mirror plane. Such analysis focuses on the instantaneous state of the crystal (dislocations lines, twins boundaries, etc), without indicating the history of the deformation.Frobenius norm of the local Green-Lagrange tensor, on the other hand, provides a measure of the difference (in terms of deformation) between the initial state (defect-free) and the actual microstructure. A dislocation loop is highlighted through the formation of a penny-shaped disregistry, representing the deformation induced by the dislocation motion. Local Lagrangian deformation gradient tensor F , when described as a composition of rotation R and stretch Û tensors (polar decomposition F = R • Û ) provides a means to distinguish zones that have undergone a reconstructive (irreversible and dissipative) deformation from zones that have only rotated (i.e. reversible deformation with no dissipation).

Figure 3 . 5 :

 35 Figure 3.5: Block-sliding configuration for the computation of E gsf on a glide plane of normal n with f the displacement jump vector on the surface S.

Figure 3 . 6 :

 36 Figure 3.6: (left) Classical representation of an edge dislocation with Burgers vector b in a plane of normal n. (right) Modelisation of the dislocation with the PNG method as a deformation field ε in an area of height L. Dislocation core spreading is characterized through the evolution of η.

3. 5 . 1

 51 Dynamique moléculaire classique 3.5.1.1 Premiers principes Formulations en Lagrangien et Hamiltonien La méthode de la dynamique moléculaire classique (DM) consiste à décrire l'évolution temporelle d'un ensemble de N particules (atomesmolécules) qui interagissent les unes avec les autres. Les équations de la DM sont basées sur une formulation Lagrangienne de mécanique classique. On considère un ensemble de N particules de coordonnées et vitesses généralisées {q i , qi }, qui décrivent l'évolution d'un système conservatif dont les forces dérivent d'un potentiel U . Le Lagrangien L = L({q i }, { qi }, t) satisfait les équations de Lagrange :

34 )

 34 Deux versions de ce potentiel ont été proposées parBedrov et al., une polarisable sans interaction intra-moléculaire entre atomes non liés et une non polarisable. Beaucoup de paramètres de ce champ de forces ont été ajustés via des calculs de chimie quantique. Tous les paramètres du potentiel original sont donnés dans le Tableau 3.1 et plus de détails sont disponibles dans l'article original[START_REF] Bedrov | A molecular dynamics simulation study of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene as a function of pressure and temperature[END_REF].Dans ce travail, nous avons utilisé la version non polarisable (qui ne tient pas compte de la polarité des interactions) du potentiel de Bedrov et al. avec un scaling de 20% sur les charges partielles, comme recommandé par les auteurs en guise d'ajustement empirique. Nous avons choisi d'effectuer des simulations de DM en utilisant des molécules de TATB entièrement rigides et nous justifions ce choix dans la partie suivante. Les interactions électrostatiques à longue distance ont été calculées via l'approximation Reaction Field[START_REF] Onsager | Electric moments of molecules in liquids[END_REF][START_REF] Barker | Monte carlo studies of the dielectric properties of water-like models[END_REF][START_REF] Watts | Monte carlo studies of liquid water[END_REF] avec un rayon de coupure de 13 Å comme pour toutes les autres interactions. Un thermostat de Langevin avec une constante de couplage de 1, 0 ps a été utilisé tandis que le pas de temps des simulations NVT et NVE a été fixé à 1, 0 fs. Enfin, les équations du mouvement sont discrétisées via l'algorithme de Verlet-Vitesse[START_REF] Verlet | Computer experiments on classical fluids. ii. equilibrium correlation functions[END_REF].

  tenseur d'élasticité C a été obtenu en température et pression en approchant un panel de courbes d'énergie élastique avec un polynôme du second degré. Ces courbes d'énergie élastique ont été calculées en appliquant des déformations successives à une boîte de simulation 3Dpériodique de taille 20a × 20b × 20c contenant 16000 molécules. Puisque le TATB est un cristal triclinique, 21 constantes élastiques ont été calculées en appliquant 21 types de déformation différents. Résultats à 300 K et 0 GPa. Le tenseur C complet est renseigné dans le Tableau 3.3. L'anisotropie du TATB apparaît clairement dans la symétrie du tenseur qui est en accord avec la structure cristalline du matériau et le type d'interactions inter-moléculaires qui intervient de manière très directionnelle. Ainsi, les coefficients de cisaillement C 44 et C 55 visant à faire glisser les feuillets de molécules les uns sur les autres sont très bas par rapport à C 66 qui traduit un cisaillement dans les feuillets de molécules. De plus on s'aperc ¸oit que le TATB est plus compressible selon la normale au plan de base que dans la direction des plans de molécules (C 33 < C 11 , C 22 ). Ces constantes élastiques sont comparables aux résultats précédemment publiés

  35) o ù α et β sont des fonctions dépendantes du temps, ⊗ est le produit tensoriel (a ⊗ b = a i b j ), (m, n) des vecteurs utilisés pour la construction de F (voir Figure 3.4) et I est le tenseur identité. Des conditions de déformation à volume et taux de déformation constants permettent d'aboutir aux expressions des fonctions α et β. Pour chaque simulation de DM, une trajectoire de 100 ps dans l'ensemble NVT est effectuée dans un premier temps. La configuration finale de cette simulation est ensuite utilisée comme point de départ pour les simulations avec déformation (t = 0). La boîte de simulation est alors déformée de manière dynamique en appliquant F aux vecteurs de périodicité (a, b, c) uniquement (voir Figure 3.4). Le champ de vitesse appliqué à l'instant initial est parfaitement cohérent avec le tenseur gradient des vitesses L déduit de F (i.e., v(t = 0) = L • x(t = 0)). Aucun rescaling des positions atomiques n'est effectué dans cette méthode, ce qui permet un mouvement cohérent des particules avec l'évolution dynamique des conditions périodiques de la boîte de simulation, et empêche la génération d'ondes acoustiques dans le matériau.

3. 5 . 3

 53 Modèle de Peierls-Nabarro-Galerkin pour les dislocations Le modèle de Peierls-Nabarro-Galerkin (PNG) est une méthode utilisée pour calculer des dislocations au coeur très complexe, qui se séparent en plusieurs dislocations partielles et pouvant être soumises à un chargement dynamique [27, 75]. Ces calculs sont en général effectués à température nulle ce qui est cohérent avec la définition des contraintes de Peierls pour les dislocations (la contrainte de Peierls est la contrainte nécessaire pour faire bouger la dislocation dans le réseau cristallin à 0 K). En revanche, la fréquence de vibration d'une dislocation possédant un coeur étalé (voir l'étude pour le TATB, ci-dessous) est plusieurs ordres de grandeur plus faible que la fréquence de vibration atomique. Une séparation en temps peut alors être introduite entre la structure de coeur de dislocation et les positions des molécules à l'équilibre. Ceci ouvre la possibilité de définir une γ-surface (l'énergie de faute d'empilement E gsf (f ) qui représente le co ût d'un saut de déplacement f dans un plan de glissement, voir Figure 3.5) en température, pouvant être utilisée dans un calcul PNG de structure de coeur de dislocation. Le modèle de PNG couple une méthode de Peierls-Nabarro (PN) généralisée avec une technique de type element-free Galerkin, comparable à la méthode des éléments finis. Cela consiste à minimiser une fonctionnelle d'énergie totale séparée en une contribution purement élastique et un terme d'énergie de glissement calculé à partir d'une γ-surface, afin d'aboutir à la structure de coeur de dislocation. Dans le modèle de PNG, le terme inélastique provient de la γ-surface de laquelle la partie élastique a été retirée. Deux champs de déplacement sont introduits : un champ 3D u et un champ 2D η contenu dans le plan de la γ-surface. Une dislocation telle qu'interprétée par le modèle PNG est représentée en Figure 3.6. Une solution pour la structure de coeur de dislocation est finalement obtenue par la minimisation de l'énergie E relativement à u et η :

  Figure 4.1 is presented a view of TATB single crystal along a. Due to the stacking of TATB molecules, two planes (denoted as P1 and P2) play the role of the (001) basal plane. Interplanar interactions between molecules are of van der Waals type while interactions between molecules belonging to the same layer are governed by strong hydrogen bonding.

[ 62 ] 2 Figure 4 . 1 :

 62241 Figure 4.1: View of the TATB single crystal along a. Glide planes of interest are denoted P1 and P2. Carbon atoms are colored in grey, nitrogen in blue, oxygen in red and hydrogen in white. The prime ( ′ ) symbol on the vector labels means that these vectors do not lie in the plane of the figure.

  MD simulations were used to provide generalized stacking fault energies and elastic constants, as input parameters of the PNG method. The fourthorder elastic tensor and the γ-surfaces have been computed for several (T , P ) conditions in order to predict the structure of edge and screw dislocations along the two Burgers vectors corresponding to the lattice vectors a and b.The chapter is organized as follows. Section 4.2 contains results from MD simulations of TATB under various conditions. The methodology proposed to obtain the γ-surfaces and the fourth-order elastic tensors at finite temperature are presented. Section 4.3 is then dedicated to mesoscopic simulations using the PNG model to study the dislocation core structure in TATB. Structures for edge and screw dislocations along the two directions within the two glide planes taken in consideration are investigated.

1 .Figure 4 . 2 :Figure 4 . 3

 14243 Figure 4.2: General procedure to compute the γ-surfaces of (001) TATB basal plane. (left) First step: equilibration in the NVT ensemble. (right) Second step: splitting of the 3D periodic simulation cell in two parts and application of an instantaneous glide of vector f .

Figure 4 . 3 :

 43 Figure 4.3: γ-surfaces for TATB at 0 K and 0 GPa for P 1 (a) and P 2 (b) planes. 300 K P 1 plane γ-surfaces at 0 GPa (c) and 10 GPa (d). MEP along [100] and [010] directions obtained with the string method are represented on the bottom left image.

Figure 4 . 4 :

 44 Figure 4.4: Energies along MEP normalized by shear elastic coefficients along [100] (a) and [010] (b) for the studied range of pressures. Curves have been shifted to highlight their similar evolution.

  norm of the gradient of η, |∇η(x)| (in Figure 4.5), which represent the density of Burgers vector. Thus, each peak of |∇η(x)| represents a partial dislocation, the width of which is characterized by the Full Width at Half Maximun (FWHM). In all simulations, dislocations appear to split into several partials, with FWHM around 10 Burgers vectors.

Figure 4 . 5 . 4 . 3 . 2 [Figure 4 . 5 :

 4543245 Figure 4.5: Frobenius norm of | ∇η | through the dislocation cores for edge dislocation with [100] (a) and [010] (b) Burgers vector and for screw dislocations with [100] (c) and [010] (d) Burgers vector.

Figure 4 . 6 :

 46 Figure 4.6: Maximum spreading length ∆ SL in number of Burgers vectors for dislocation core structures at 300 K. (Continuous lines) Edge dislocations. (Dashed lines) Screw dislocations.

  and the decrease of the interplanar distance results in an increase of the γ-surfaces energy barriers. Dislocations are predicted to have complex core structures within P 1 and P 2 planes and to spread over approximately a hundred of Burgers vectors. A splitting into up to four partials (with a width at half-maximum around 10 lattices), regardless of the plane, is observed for both[START_REF] Stevens | Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder x-ray diffraction[END_REF] and [010] Burgers vectors. Screw dislocations are approximatively 25% thinner than edge ones. Thus, dislocation friction on the lattice is unlikely to play a dominant role for the viscoplastic behavior of TATB material for the range of pressures and temperatures explored.The estimation of the maximum von Mises stress along the MEPs at ambient conditions is around 250 MPa which defines an upper bound for the shear stresses within (001) planes.

Figure 4 . 1 . 4 . 5 . 1 4 . 5 . 1 . 1 2 .

 4145145112 En raison de l'empilement des molécules, deux plans (notés P1 et P2) de base (001) peuvent être distingués. Les faibles interactions de van der Waals qui régissent les interactions inter-plans constrastent avec les fortes liaisons hydrogène présentes dans les feuillets, indiquant que le glissement de dislocations dans les plans P1 et P2 doit être très favorable. Les dislocations, connues pour être dominantes dans le comportement plastique des alliages métalliques et minéraux, n'ont jamais été observées dans le TATB. Les faibles énergies de glissement pourraient mener à des dislocations au coeur très étalé et donc la nécessité de grandes boîtes de simulation en dynamique moléculaire. Afin de contourner ce problème, le modèle de Peierls-Nabarro Galerkin (PNG) présenté en partie 3.3 est utilisé. Le but de ce chapitre est la détermination des structures de coeur de dislocations dans le plan basal du TATB en température et en pression. Des simulations de DM sont utilisées pour calculer les γ-surfaces qui sont utilisées avec les constantes élastiques déterminées au chapitre précédent dans le modèle PNG. Une des nouveautés de ce travail est le calcul en pression et température. Ces deux propriétés ont été calculées pour différentes pressions entre 0 et 10 GPa à 0 et 300 K afin de caractériser les dislocations coin et vis de vecteurs de Burgers [100] et [010]. Propriétés du glissement dans les plans P1 et P2 Calcul des γ-surfaces Une γ-surface est une énergie potentielle calculée en appliquant un glissement f dans un plan de glissement d'un cristal parfait. On déplace ainsi de f le demi cristal supérieur par rapport au demi cristal inférieur, en n'autorisant que des mouvements selon la normale au plan de glissement dans les deux demi cristaux. Ces calculs sont généralement effectués à 0 K avec une relaxation des atomes perpendiculairement au plan de glissement. Une définition des γ-surfaces en température a été proposée dans le but de prendre en considération l'effet de la température sur les dislocations dans le TATB. Une telle mesure en température nécessite en général une moyenne temporelle. En revanche dans le cas d'une γ-surface, un équilibrage en température n'est pas possible car la faute d'empilement provoquerait une forte modification de l'arrangement atomique local. Pour contourner ce problème, nous avons équilibré un échantillon à 300 K avant d'appliquer la déformation et de récupérer l'énergie totale moyennée sur toutes les particules du système. Par ergodicité, cette ensemble de configurations équivaut à un grand système. Un cristal de taille 20a × 20b × 20c est considéré par la suite. Une base non orthonormale (a, b, c) définit la maille cristalline. Une seconde base orthonormale (x, y, z) est définie avec a parallèle à x, b dans le plan (x, y) et c orienté dans le demi espace supérieur. Deux γ-surfaces doivent être calculées pour étudier les deux plans P 1 et P Après une trajectoire dans l'ensemble NVT de 100 ps, le cristal est séparé en deux dans le plan P 1 (respectivement P 2) et un glissement instantanné est appliqué dans le plan (x, y) (voir Figure 4.2). Afin d'éviter le recouvrement d'images périodiques dont le déplacement est différent et de calculer une énergie cohérente, la périodicité du cristal est supprimée dans la direction z. Une telle γ-surface n'est pas parfaitement périodique puisqu'après un déplacement d'un vecteur de maille, l'énergie n'est pas la même qu'à déformation nulle à cause des différentes orientations et positions de chaque molécule. Afin d'atténuer cet écart, nous avons en fait choisi de calculer 16 γ-surfaces puis de les moyenner sur l'intervalle f ∈ [-2a, +2a], [-2b, +2b]). Ceci réduit la différence en énergie de deux ordres de grandeur, à moins de 0, 5 %. Les γ-surfaces à 0 K et 300 K et pour 5 pressions entre 0 et 10 GPa ont été calculées. Les γ-surfaces à 0 K pour les plans P 1 et P 2 sont données en

Figure 4 . 3 ; 4 . 5 . 1 . 2 4 . 5 . 2 4 . 5 . 2 . 1 4 . 5 . 2 . 2 1 , 5 . 4 . 5 . 3 Conclusion

 4345124524521452215453 Figure 4.3 ; elles sont similaires à celles calculées par Mathew et Sewell [73]. Les résultats à 300 K montrent que la température adoucit le paysage énergétique, sans effet majeur sur les énergies de faute d'empilement. 4.5.1.2 Équilibre entre énergie élastique et inélastique La pression ne modifie pas la forme des γ-surfaces, ceci étant d û à la grande stabilité de la structure triclinique du TATB et à la conservation de son anisotropie élastique. Cette propriété est vérifiée en calculant les chemins minimum d'énergie (MEPs) pour des glissements [100] et [010] en utilisant la méthode de la corde [113] pour les différentes pressions (voir Figure 4.4). Ces courbes d'énergie ont été normalisées par le module de cisaillement moyen dans le plan (B 44 + B 55 )/2, estimé à partir des coefficients de Birch, afin d'étudier l'équilibre entre énergie élastique et inélastique. On s'aperc ¸oît que les courbes suivent la même évolution, ce qui montre que la structure de coeur des dislocations ne devrait pas être très impactée par la pression. La présence de plusieurs minima le long des MEPs indique que la dislocation aura tendance à se séparer en plusieurs partielles. Nos résultats concernant les γ-surfaces et les constantes élastiques semblent être cohérents avec ceux de la littérature, ce qui nous permet de les utiliser dans le modèle de PNG.4.5.2 Structure de coeur des dislocationsUne boîte de simulation équivalente à 1024 × 256 mailles de TATB est considérée. Une dislocation parfaite est introduite dans le cristal avant de minimiser l'énergie totale du système.4.5.2.1 Dislocations[100](001) et [010](001) La structure du coeur des dislocations est représentée en utilisant la norme de Frobenius du gradient du déplacement |∇η(x)| (voir Figure 4.5), qui représente l'intensité du vecteur de Burgers. Chaque pic correspond à une dislocation partielle dont l'étalement est caractérisé par la largeur à mi-hauteur (FWHM) du pic. Dans toutes les simulations, les dislocations se séparent en plusieurs partielles avec une valeur de la FWHM de 10 vecteurs de Burgers environ. Le comportement des dislocations dans les

Figure 5 . 1 : 1 . 5 . 2 . 2 Figure 5 . 2 : 1 to 3 .

 5115225213 Figure 5.1: Stress/strain curve for a compression along direction m with θ = 30 • and φ = 0 • at a strain rate of 1 × 10 8 s -1 .

Figure 5 . 3 .

 53 Figure 5.3. This nucleation of dislocations in transverse planes always comes with an unusual local dilatancy process, that will be discussed in the next paragraph.

Figure 5 . 3 : 5 . 3 . 1

 53531 Figure 5.3: TATB single crystal non-basal gliding in the (0 11) plane for a deformation that consists in a pure shear along this plane. The prime ( ′ ) symbol on the vector labels means that these vectors do not lie in the plane of the figure. All-atom display with identification of basal (P 1,P 2) and transverse planes (0 11) of interest.

Figure 5 . 4 :

 54 Figure 5.4: Review of potential slip vectors in transverse planes of TATB single crystal. (a) View along a. (b) View along b. Brown and orange shapes represent TATB molecules. Orange molecules are in the figure plane whereas brown ones are located behind it. Burgers vectors are represented as green continuous arrows and partial slip vectors as black dashed arrows. The prime ( ′ ) symbol on the vector labels means that these vectors do not lie in the plane of the figure.

  Figure 5.5a, red annotation), which comes with a local strain that transforms the hexagonal stacking (yellow annotation in Figure 5.5b) into a rectangular-like arrangement. For the plane with the lower dilatancy, a strain is also visible, although less stronger (blue annotation in Figure 5.5c).

Figure 5 . 5 :

 55 Figure 5.5: Transverse dislocation process. (a) 100 view of three successive planes during the midway configuration of a [011](011) dislocation that propagates into a TATB single crystal. Blue and red shapes encircle stacking fault configurations SF 1 and SF 2 respectively. (b) Molecular arrangement in a defect-free TATB basal plane. (c) and (d) In-plane molecular stacking for the two successive stacking fault configurations in TATB basal plane SF 1 and SF 2, respectively. Carbon is cyan, nitrogen is blue, oxygen is red and hydrogen is white.

Table 5 . 1 :

 51 Activated Slip Systems in TATB Single Crystal along the (001) plane. We address this point in this chapter, starting by a definition of the geometrical transformation for the twinning in section 5.4.1. Since both mechanisms involve pure shear of the basal plane, the energy landscape for this deformation may have an influence on the resulting microstructure. In section 5.4.2, the computation of this energy landscape that details the energetics during any simple shear in TATB basal plane is proposed. Finally, a study of large TATB single crystal plates (100a × 100b × 10c MD simulation cells) is undertaken in section 5.4.3

2 Figure 5 . 6 :

 256 Figure 5.6: Buckling mechanism of TATB single crystal for a deformation that mainly contains a compressive component along the [100] direction (θ = 0 • , φ = 0 • ). The prime ( ′ ) symbol on the vector labels means that these vectors do not lie in the plane of the figure. All-atom display with identification of basal (P 1,P 2) planes.

Figure 5 . 7 :

 57 Figure 5.7: Detailed twinning process for TATB single crystal. (a) Three unstrained TATB triclinic cells. (b) TATB cells after inversion of the two out-of-plane angles α and β. The unit cell shifted about 1 2 c (highlighted in red) is the original one after a rotation of 180 • around the z-axis. Only the process of the α angle inversion is shown here, which is similar for the β angle. Carbon is cyan, nitrogen is blue, oxygen is red and hydrogen is white.

Figure 5 . 8 :

 58 Figure 5.8: Homogeneous shear energy surfaces for TATB basal plane. (a) Shear energy as a function of the right Cauchy-Green tensor C components: E shear = g(C 13 , C 23 ). (b) Shear energy as a function of out-of-plane angles: E shear = g(α, β).

EFigure 5 . 9 :Figure 5 . 9 .

 5959 Figure 5.9: (left) Energies along the two minimum energy paths (MEPs) linking the stable and twin configurations of TATB. (right) MEPs reported on the twinning energy landscape.

Compression with θ = 30 • 1 × 10 10 s - 1 2 × 10 9 s - 1 1 × 10 9 s - 1 5 × 10 8 s - 1 1 × 10 8 Figure 5 . 10 :

 30111118510 Figure 5.10: Effect of strain rate on the nucleation stresses. (left) Strainrate curves for a compression along direction m with θ = 30 • . (right) Different values of the nucleation von Mises stress at nucleation as a function of strain rate, for different directions of compression.

Figure 5 .

 5 10 are presented the stressstrain curves for a loading along m with θ = 30 • and strain rates in [1 × 10 8 s -1 , 1 × 10 10 s -1 ] showing a decrease of the stress to nucleation for decreasing strain rate. No effect of the compression angle θ was observed for the nucleation stress.

Figure 5 . 11 :

 511 Figure 5.11: Strain field evolution with respect to the strain rate. The final microstructure for 4 different molecular dynamics simulations are shown. Each snapshot corresponds to a compression along m with θ = 30 • . Color code corresponds to the value of ||E||.

Figure 5 . 12 :

 512 Figure 5.12: Microstructure evolution with respect to the compression direction. (rows) Snapshots for three different values of the angle θ. (columns) Snapshots for three different values of the strain along the prescribed deformation path simulation. Color code corresponds to the value of E.

E 23 -Figure 5 . 13 :

 23513 Figure 5.13: Link between the homogeneous shear energy surface and the deformation signatures obtained during directional deformations. The gray color bar represents the shear energy (similar to Figure 5.8) while the filled circles represents the deformation signatures, each color corresponding to a different direction of compression.

EFigure 5 . 14 :

 514 Figure 5.14: Energies of the twinning energy landscape along different directions passing by the α-TATB phase minimum.

Figure 5 . 15 : 5 . 5 Pressure, size and strain rate effects 5 . 5 . 1

 51555551 Figure 5.15: Evolution with pressure of the nucleation von Mises stress σ v (φ) (θ = 90 • ). Gray triangles correspond to the simulation data, curves are interpolation to guide the eye.

Figure 5 .

 5 Figure 5.16: (a) All-atom representation of six strain-free TATB layers. The buckling mechanism can be decomposed into a stretch deformation (b) followed by a rigid body rotation (c).

Figure 5 . 16 . 7 )

 5167 The local entire deformation is then governed by a deformation gradient tensor F = R • U with R a rotation matrix of angle θ and Read et al. define F = I + ∇v where v = (v x , v y , v z ) is the displacement field of the layered material. Replacing F by its analytical expression leads to a system that can be resolved to express θ, ε 22 , ε 23 and ε 33 as functions of the displacement gradient ∇v:

Figure 5 . 17 :

 517 Figure 5.17: Colormap of the function f 2D (ε yy , εzz ). Negative values are orange, positive are black. The black thick line represents f 2D (ε yy , εzz ) = 0. The black circle represents the function at zero deformation, the star at (ε yy , εzz ) = (-0.043, 0.0) and the square at (ε yy , ε zz ) = (-0.022, 0.05).

2 . 1 , 5 . 8 . 1 . 1 •

 215811 construite pour imposer une déformation suivant un chemin est utilisée pour déclencher la nucléation de défauts dans le monocristal de TATB. Des courbes contrainte-déformation pour des chargements de cisaillement et de compression uniaxiale sont obtenues pour plusieurs directions. De plus, la mesure locale de la déformation présentée en partie 3.2.2 (détaillée dans l'annexe A.3) permet d'analyser finement les défauts nucléés. 5.8.1 Surface d'écoulement plastique pour le monocristal de TATB Des simulations de DM sont effectuées avec un système de taille 10a × 10b × 14c. La structure cristalline pour différentes températures et pressions a été obtenue par des simulations DM de Parrinello-Rahman. Le système est dans un premier temps équilibré dans l'ensemble NVT, puis une déformation de cisaillement pur lui est appliquée en suivant le formalisme présenté en partie 3.2.1. Le plan défini par les vecteurs m et n est choisi perpendiculaire au plan de base du TATB. Deux conditions sont appliquées pour ces simulations : les déformations préservent le volume et sont effectuées à taux de déformation constant. Le tenseur gradient de la déformation dépendant du temps t est défini dans l'équation 5.1. Les vecteurs m et n sont perpendiculaires et représentent une direction de compression et de dilatation. Les coefficients devant chaque produit tensoriel sont établis, afin que la déformation soit exactement égale à un cisaillement pur dans la direction m + n dans un plan de normale nm. Résultats à (300 K, 0 GPa) Les mécanismes de déformation du TATB ont été caractérisés à 300 K et 0 GPa. 84 couples d'angles θ et φ ont été choisis afin de solliciter le TATB en cisaillement pur dans toutes les directions de l'espace. La déformation équivalente finale de chaque simulation est égale à 10%. Pour chaque trajectoire, une chute de la contrainte globale du système a été observée et interprétée comme la nucléation du premier défaut. La contrainte de nucléation σ v (θ, φ) est définie par la contrainte de von Mises juste avant nucléation, fortement dépendante de la symétrie du TATB. Une surface 3D de σ v (θ, φ) est présentée en Figure 5.2 dans laquelle 3 zones distinctes apparaissent : La zone 1 correspond à φ > 45 • . Le glissement dans le plan basal est très peu activé et la déformation est accommodée par la nucléation de dislocations dans des plans non basaux. La nucléation de tels défauts implique une rupture des liaisons hydrogène dans les feuillets des molécules ainsi qu'une forte dilatance locale du coeur de la dislocation. • La zone 2 correspond à la plus faible valeur de la contrainte de nucléation. La déformation est effectuée pour φ = 45 • et correspond à un cisaillement pur dans le plan de base du TATB. Ce glissement est activé pour une contrainte inférieure à la précision du calcul de la contrainte, en raison des faibles interactions de van der Waals entre les feuillets de molécules. • La zone 3 quant à elle concerne les déformations pour lesquelles φ < 45 • et o ù la composante de compression se trouve dans la direction des plans de molécules. Ces déformations donnent lieu à un mécanisme de flambage, provoquant une forte désorientation des feuillets de leur plan initial.

Figure 5 .

 5 Figure 5.4 contient deux vues des plans (100) and (010), orthogonaux à b × c et c × a. Les flèches continues vertes représentent des vecteurs de Burgers tandis que les flèches noires servent à identifier des vecteurs de

Figure 5 .

 5 6). Une telle déformation peut être vue comme une composition d'un cisaillement pur et d'une rotation. Pour le TATB, lors de l'activation de ce mécanisme, le plan (001) subit un fort cisaillement. Nous avons dans un premier temps étudié les conséquences géométriques d'un cisaillement homogène du plan basal du TATB. La symétrie de type hexagonale du plan basal et l'empilement des molécules permettent d'obtenir un maclage parfait du monocristal lorsque celui-ci subit un cisaillement pur dans le plan (001) qui inverse les deux angles hors plan α et β par rapport à 90 • (α ′ = 180 -α = 91, 66 • , β ′ = 180 -β = 76, 66 • ). Cette déformation change l'orientation des groupes amino et nitro au sein de la maille déformée ainsi que les positions des centres de masse des molécules. Si on considère la maille décalée de 1 2 [001] (voir Figure 5.7), on restaure la maille de TATB initiale, qui a été tournée de 180 • . Afin de vérifier cette hypothèse de maclage parfait du TATB, nous avons appliqué des chemins de déformation pour un cisaillement homogène dans le plan de base afin de calculer une surface d'énergie de maclage. Le tenseur gradient de la transformation utilisé pour cette déformation n'a que deux composantes non triviales F 13 et F 23 . En prenant le tenseur de Cauchy-Green correspondant (C = F T • F ) afin d'éliminer les rotations pures, une énergie potentielle d'excès est calculée en fonction de C 13 et C 23 . Ce paysage énergétique est représenté en Figure 5.8a, sur lequel le minimum d'énergie correspondant à l'état stable du TATB est noté α-TATB. Dans le voisinage de ce minimum stable, on distingue un second minimum, qui correspond exactement au maclage de la maille cristalline, noté α T -TATB. On peut voir en Figure 5.7b que ce minimum correspond bien à l'inversion des angles α et β. En utilisant la méthode de la corde comme pour les γ-surfaces, on peut calculer la contrainte moyenne le long des chemins de minimum d'énergie. On montre ainsi que la contrainte maximale de cisaillement nécessaire pour provoquer le maclage n'est que de l'ordre de 150 MPa. Afin d'étudier le mécanisme de flambage et de le relier au cisaillement nécessaire pour le maclage, nous avons étudié le comportement du TATB en compression simple dans la direction des feuillets de molécules. La taille de la boîte de simulation utilisée est de 100a × 100b × 10c. Deux types de simulations sont effectuées : le premier consiste à comprimer de manière instantannée le matériau afin d'observer la dépendance du comportement par rapport à la déformation appliquée et à la direction de compression. Le second comprend une première trajectoire NVT puis une phase de déformation à taux de déformation constant jusqu'à 5% de déformation ce qui permet d'étudier la dépendance de la microstructure par rapport à la vitesse de déformation. La dépendance à la vitesse de déformation a été étudiée pour des trajectoires de déformation entre 50 ps et 500 ps. En Figure 5.10 sont présentées les courbes contrainte-déformation pour une compression selon m avec θ = 30 • et une vitesse de déformation couvrant l'intervalle [1 × 10 8 s -1 , 1 × 10 10 s -1 ]. On observe globalement des contraintes de nucléation plus faibles lorsque l'on baisse la vitesse de déformation. De plus, aucun effet sur la valeur de la contrainte de nucléation n'a été observé en changeant la direction de compression. En Figure 5.11 sont présentées des images de la microstructure o ù E est tracé par molécule pour les microstructures finales obtenues à différentes vitesses de déformation. E = 1 2 (F T F -I) est le tenseur de Green-Lagrange calculé à partir du tenseur gradient de la transformation locale. On observe de manière très claire que l'augmentation de la vitesse de déformation provoque un nombre plus élevé de sites de nucléation et une microstructure plus perturbée. De plus, la taille des bandes de déformation (appelées "chevrons") présente également une dépendance, avec des chevrons plus larges à faible vitesse de sollicitation. Cependant, la taille finale de ces microstructures dépend probablement de la taille du système. En effet, une fois que la totalité de l'interface entre deux chevrons est parallèle à une direction de périodicité, la microstructure a atteint l'équilibre et n'est plus modifiée. Plusieurs simulations à la même vitesse de déformation de 2×10 8 s -1 ont été effectuées afin d'étudier la dépendance de la microstructure à la direction de compression. La boîte de simulation a ainsi été comprimée selon m = (cos θ, sin θ, 0) avec θ = 30, 75 et 90 • . Les chevrons sont systématiquement obtenus avec des interfaces perpendiculaires à la direction de compression comme montré en Figure 5.12. Les simula-tions de DM pour les trois différentes orientations de compression sont représentées pour trois valeurs intermédiaires de la déformation longitudinale. Ainsi, en partant d'un cristal parfait, les premières désorientations apparaîssent pour environ 3% de déformation, avec une nucléation homogène de petits chevrons. Au fur et à mesure que la déformation progresse, un processus de coalescence hiérarchique opère avec le regroupement de plusieurs chevrons en bandes de déformation plus larges. Le seuil de nucléation du flambage a été estimé via des simulations à différentes valeurs de déformation longitudinale, suivies d'une trajectoire NVT afin de voir (éventuellement) émerger une instabilité de flambage. Ce seuil de nucléation est noté T B et vaut T B = 2, 75% ± 0, 25 pour toutes les directions de compression. Puisque le flambage induit un fort cisaillement dans le plan basal, les seules composantes non triviales du tenseur de Cauchy-Green calculé localement sont C 13 and C 23 . L'analyse du couple (C 13 , C 23 ) est importante si l'on veut relier le flambage à la surface d'énergie de maclage calculée précédemment qui est tracée en fonction de ces même composantes du tenseur de Cauchy-Green. Ainsi le couple (C 13 , C 23 ) est calculé pour chaque molécule afin de définir une carte de densité d(C 13 , C 23 ), telle que le nombre total de molécules pour l'intervalle de déformation (C 13 ± ∆ 2 , C 23 ± ∆ 2 ) est donné par d(C 13 , C 23 )×∆ 2 . Cette carte de densité est ensuite analysée afin d'en obtenir les positions des différents maximums, correspondantes aux différentes zones ayant subies un cisaillement lors des simulations. En reportant ces différents couples des composantes du tenseur de Cauchy-Green sur le paysage énergétique de maclage, on peut déterminer les éventuelles corrélations entre ce paysage et la microstructure. Pour les trois simulations discutées plus tôt, les histogrammes simplifiés sont tracés en Figure 5.13. Pour une compression avec θ = 30, la déformation est principalement en C 13 et le couple de déformation est aligné avec la direction de compression. De manière similaire, lorsque θ = 90 • , le couple est aligné avec la direction de compression mais cette fois-ci la déformation est en C 23 . Concernant le cas o ù θ = 75 • , le point le plus éloigné du minimum énergétique stable est séparé en deux et une partie se trouve dans le voisinage du minimum correspondant au maclage. Ces observations montrent que lorsque la direction de compression est proche des deux minima, la microstructure est fortement dépendante du paysage énergétique de maclage. En revanche, plus la direction de compression se rapproche du cas θ = 0 • , plus la microstructure est déterminée par un phénomène de flambage simple. Nous avons appelé ce comportement mixte "Maclage-Flambage" dans la suite. Certains couples tombent dans le voisinage des minima locaux, confirmant ainsi le rôle du paysage énergétique sous-jacent. En comparant l'excès en énergie le long des MEPs (estimé à 25, 4 MJ m -3 ), et l'énergie élastique disponible avant l'apparition du flambage, on peut estimer qu'une simple compression peut mener à l'obtention du maclage pour le monocristal de TATB. L'énergie élastique pour une déformation de

  ments distincts. Le glissement dans le plan basal est très facile à activer et est relié aux faibles énergies de glissement dans ces directions. Ceci est cohérent avec les résultats du chapitre précédent. Des dislocations transverses ont été identifiées et représentent un nouveau mécanisme de plasticité pour le monocristal de TATB. Cependant, leur comportement semble très complexe avec une forte dilatance locale du coeur et un processus de nucléation impliquant des ruptures de liaisons hydrogène dans les feuillets de molécules. Finalement, lorsque la déformation comprend une composante de compression dans les feuillets de molécules, un phénomène de flambage apparaît systématiquement à partir de 3% de déformation. La microstructure résultante consiste en une succession de chevrons de déformation opposée avec des interfaces perpendiculaires à la direction de sollicitation. Le calcul d'un paysage énergétique de cisaillement homogène dans le plan basal a également été effectué. Il a révélé un second minimum énergétique correspondant à un mécanisme de maclage parfait de la maille cristalline. La transition entre le flambage et le maclage a été démontrée comme fortement dépendante de ces énergies ainsi que de la direction de compression subie par le matériau. De plus, nous avons montré que le flambage est une instabilité purement élastique, à l'aide d'une analyse théorique basée sur des travaux antérieurs. Tous ces éléments sont des ingrédients nécessaires à la construction d'une loi de comportement aux échelles supérieures qui viserait à reproduire idéalement la plasticité via les dislocations, le maclage parfait et le mécanisme de flambage. The results of the precedent chapters have led to a much better understanding of TATB single crystal behavior. Various new deformation mechanisms have been revealed, depending on the direction of loading:

44 Figure 6 . 2 :

 4462 Figure 6.2: Different EOS computed through MD simulations. (triangle marks) MD data from the hydrostatic pressure conditions. (circle marks) MD data from the spherical deformation conditions. (red line) Vinet fit on the hydrostatic pressure data. (orange line) Vinet fit on the spherical deformation data. (grey line) Vinet fit on the spherical deformation after the pseudo phase-transition. The insert contain the parameters of the different Vinet fits.

  .3a is shown a snapshot of the system with three different color codings corresponding to three deformation estimations at each molecular center of mass. The first snapshot is the norm of the local Green-Lagrange strain tensor E = 1 2 (F T F -I). Two kind of deformed zones appear. The ones with interfaces perpendicular to x and the others with interfaces perpendicular to y. They undergo different strain level. The second view is colored with respect to the local value of C 13 that stands for pure shear in the x direction within the layers. Each vertical interface separates zones of deformation with same amplitude but opposite signs which is clearly reflected in 6.3c where the slice (A-A') has been reported. The zones with horizontal interfaces (visible on the first view) can however no longer be seen. Conversely, the third view is colored with respect to the local value of C 23 that stands for pure shear in the y direction. This time, the horizontal interfaces are clearly visible and separate zones with deformation of opposite signs but different amplitude. The slice (B-B') is reported in Figure 6.3 where the disequilibrium between successive chevrons plainly emerges.

Figure 6 . 3 : 13 Figure 6 . 4 :

 631364 Figure 6.3: (top row) Deformation snapshots of the TATB sample during MD simulations of size 100a × 100b × 10c. (left) Color code corresponds to the value of the Green-Lagrange equivalent strain ||E||. (middle) C 13 component of the Cauchy-Green strain tensor. (right) C 23 component of the Cauchy-Green strain tensor. (bottom row) View of the slices (A-A ′ ) corresponding to pure buckling (left) and (B -B ′ ) corresponding to pure twinning.

  code was used in a 3D Lagrangian framework, with an explicit time integration. Different material parameters were used for the simulation of the TATB single crystal behavior: density ρ 0 based on the single crystal equilibrium volume obtained via MD Parrinello-Rahman simulations, fourth-order elastic tensor at (300 K, 0 GPa) C 0 presented in chapter 3 and a Vinet EOS fitted on the data obtained in MD through the application of a spherical deformation with parameters K SD = 17.87 GPa, K ′ SD = 12.75 and an initial volume corresponding to density ρ 0 .

2 .

 2 Then, from t = 0.1t f to t = t f the volume is kept constant, allowing the microstructure to relax.Remark: This relaxation procedure is similar to the conditions of the MD simulation. Furthermore, since the MD simulation is equilibrated for t f = 800 ps, with a length of L M D ≈ 0.1 µm, we can compute the time needed for an acoustic wave to run through the MD sample, i.e. t c (M D) = L M D /c L ≈ 33 ps. Thus, during t f = 800 ps, the acoustic waves run N aw = t f /t c (M D) ≈ 25 times through the sample in the MD simulation. This quantitative analysis shows that the microstructure presented for the MD simulation corresponds in fact to the microstructure obtained in the Lagrangian framework at t = 0.25t f = 25t c .1

Figure 6 . 6 :

 66 Figure 6.6: TATB single crystal elastic buckling simulation. (top) Stressstrain histograms of the EB simulation. Colored points are the histograms at different instants and grey empty squares the average values of the EB simulation. Dark red and orange lines represent the two Vinet fits from MD simulations. (bottom left) Single crystal microstructure from the EB simulation. (bottom right) Single crystal microstructure from the MD simulation.

Figure 6 . 7 :

 67 Figure 6.7: Polycrystalline TATB elastic buckling simulation. (top) Stress-strain histograms of the EB simulation. Colored points are the histograms at different instants and grey empty squares the average values of the EB simulation. Dark red and orange lines represent the two Vinet fits from MD simulations. (bottom) Polycrystal microstructure from the EB simulation.

  Les résultats du chapitre précédent ont permis une meilleure compréhension du comportement du TATB monocristallin. Plusieurs mécanismes de déformation ont été révélés : les simulations de dynamique moléculaire prédisent une nucléation de dislocations hors du plan basal ainsi qu'un cisaillement pur menant à un mécanisme de maclage parfait. De plus, un mécanisme de maclage-flambage apparaît lorsque qu'une faible compression uniaxiale est appliquée aux feuillets de molécules. La microstructure résultante dépend fortement du paysage énergétique lié au maclage. Ce chapitre présente des résultats nouveaux sur le phénomène de flambage à l'aide de simulations à l'échelle mésoscopique. Lorsqu'il s'agit de construire des lois de comportement prédictives pour le comportement sous choc des matériaux, les propriétés physiques doivent être évaluées dans un certain intervalle de pressions. Parmi les propriétés d'intérêt, les relations volume-pression sont modélisées en utilisant des équations d'état. Ces équations peuvent être mesurées expérimentalement ou calculées via des simulations numériques puis approchées par des fonctions empiriques (équations d'état de Birch-Murnaghan[START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF][START_REF] Murnaghan | The compressibility of media under extreme pressures[END_REF][START_REF] Birch | Finite elastic strain of cubic crystals[END_REF] et de Vinet[START_REF] Vinet | A universal equation of state for solids[END_REF] par exemple).

Figure 6 .

 6 2 montrent qu'en présence d'une pression hydrostatique l'évolution de la relation pression-volume est compatible avec une équation d'état et les points de DM sont approchés avec une équation de Vinet (HP EOS) dont le module de raideur K HP calculé à partir du tenseur des souplesses obtenus en DM. Lorsqu'une déformation sphérique est appliquée, on observe un comportement bien différent. La relation pression-volume ne peut être approchée par une équation de Vinet que sur l'intervalle V V 0 ∈ [1.0, 0.88] avec le module de raideur K SD calculé à partir du tenseur des constantes élastiques calculé en DM. Au-delà on montre l'apparition d'une instabilité de type Maclage-Flambage similaire à ce qui a été observé dans le chapitre 5. La génération de cette instabilité est visible jusqu'à V V 0 = 0.81 puis le matériau voit ses propriétés élastiques fortement modifiées, ce qui s'apparente à une transition de phase du premier ordre. Lorsque la déformation se poursuit, on observe une évolution pression-volume à nouveau compatible avec une équation d'état. Afin de mieux comprendre la microstructure générée lors de la déformation sphérique, une simulation de DM a été réalisée avec une compression telle que V V 0 = 0.85 et une taille de système de 100a × 100b × 10c. Après 800 ps d'équilibrage dans l'ensemble NVT, la microstructure résultante est montrée en Figure 6.3. Deux types de déformations ressortent : les zones à déformation équivalente séparées par des interfaces de normale x correspondent à du pur flambage. En revanche, les zones de faible et forte déformation séparées par des interfaces de normale y correspondent à du maclage pur. Les histogrammes de la déformation correspondant à cette microstructure sont reportés sur le paysage énergétique de maclage en Figure 6.4 o ù l'on voit que les points de flambage sont centrés en zéro et en déformation C 13 uniquement alors que les points de maclage sont dans le voisinage des deux minimums avec une déformation en C 23 uniquement. 6.7.3 Résultats expérimentaux sur le TATB Un bilan des résultats expérimentaux pour le TATB a été effectué et les données sont présentées en Figure 6.5. Bien que différents types d'échantillons aient été utilisés (poudre et polycristaux) dans une cellulee à enclumes de diamant, on remarque que les équations d'état obtenues suivent toutes la tendance de l'équation d'état en pression hydrostatique alors qu'aucune de ces expériences n'a considéré un monocristal pur. La suite de ce chapitre propose une discussion et des résultats préliminaires sur l'influence de la contrainte non-hydrostatique sur le comportement à fonction de la variation de volume :

. 31 ) 6 . 7 . 4 . 2 6 . 7 . 4 . 3

 3167426743 avec C 0 le tenseur d'élasticité à l'ambiante et B 0 le module de raideur associé. Validation sur le TATB monocristallin Un code aux "element free Galerkin" (proche de la méthode des éléments finis) en Lagrangien avec intégration explicite a été utilisé prenant en entrée la densité du TATB et son tenseur d'élasticité fonction de V V 0 déduit de l'équation d'état SD EOS calculée en dynamique moléculaire pour le monocristal de TATB. Le matériau simulé comporte 64000 noeuds d'interpolation dans une boîte 3D périodique. Une simulation en deux temps est effectuée : une déformation sphérique est appliqueée jusqu'à V V 0 = 0, 85 sur un dixième du temps de calcul puis le système est relaxé à volume constant. Les histogrammes contrainte-volume sont reportés en points colorés en Figure 6.6 ainsi que leur valeur moyenne. Jusqu'à V V 0 = 0, 88 les histogrammes suivent parfaitement la courbe SD EOS confirmant que le calcul d'élasticité non linéaire est correctement effectué. Durant cet intervalle, le TATB est déformé uniformément et aucune instabilité n'est générée. Cependant, lorsque V V 0 < 0, 88, une dispersion des histogrammes est observée et les valeurs moyennes évoluent entre les courbes HP EOS (borne inférieure) et HP EOS (borne supérieure). L'instabilité élastique est donc observée pour la même compression qu'en DM puis la relation contrainte-volume évolue entre les deux équations d'état calculées en DM. Une fois que le volume final est atteint, la relaxation mène à une microstructure équilibrée en déformation. En Figure 6.6, on voit que les amplitudes de la déformation sont égales que l'on soit en DM (notée MD Simulation) ou en élasticité non linéaire (notée EB Simulation). En revanche, dans l'image provenant de la DM, 50% du volume du matériau est affecté par une déformation nulle en E 13 et correspond en fait au mécanisme de maclage. Cette caractéristique est absente dans la simulation en élasticité non linéaire car le paysage énergétique du maclage n'a pas été pris en compte. Simulations du TATB polycristallin Puisque les simulations de DM et d'élasticité non linéaire sont en bon accord qualitatif, nous avons étudié le comportement d'un polycristal de TATB en suivant une déformation sphérique afin d'étudier l'influence de la contrainte non-hydrostatique sur le polycristal. Ce dernier est créé en effectuant un pavage de Voronoï 3D périodique de l'espace comprenant 10 grains de TATB dont l'orientation est définie aléatoirement. Un seul tirage a été effectué et on pourrait envisager de considérer d'autres orientations. Une simulation en deux étapes est effectuée sur ce polycristal avec une compression de V V 0 = 0.80 avant relaxation à volume constant. Les histogrammes contrainte-volume sont reportés en Figure 6.7 avec leur valeur moyenne. On voit que les points se dispersent très rapidement à cause de l'activation de l'instabilité élastique de flambage et s'écartent fortement de la courbe SD EOS au fur et à mesure que la déformation continue mais restent situés entre les courbes HP EOS et SD EOS. Les premiers signes du flambage dans les grains apparaissent pour un taux de V V 0 = 0.88, similaire aux résultats précédents. La dispersion des histogrammes est donc due aux interactions élastiques entre grains. Du fait de la déformation macroscopique sphérique qui est appliquée au polycristal, chaque grain subit une déformation générant des contraintes non-hydrostatiques à cause de l'orientation aléatoire de ses voisins. Cependant, la déformation subie par chaque grain n'est pas forc ¸ément sphérique. En Figure 6.7, le code couleur pour la microstructure est basé sur E 13 et on voit que tous les grains subissent un flambage. 6.7.5 Conclusion Nous avons calculé une équation d'état du monocristal du TATB en pression hydrostatique via la dynamique moléculaire, en bon accord avec l'expérience et ne menant à aucun changement de phase. Une autre équation d'état en déformation sphérique a été calculée et révèle une instabilité de maclage-flambage que l'on peut assimiler à une transition de phase du premier ordre à partir de V V 0 = 0, 90. La microstructure résultante consiste en un mélange de pur flambage et de maclage. Ce résultat surprenant révèle, qu'en déformation sphérique, une pseudo transition de phase peut être observée, ce qui est impossible sous pression purement hydrostatique. L'équation d'état en déformation sphérique calculée en DM a ensuite été utilisée dans des simulations d'élasticité non linéraire dans un contexte Lagrangien. Une déformation sphérique a été appliquée à un monocristal de TATB, menant à un bon accord qualitatif entre les deux techniques de simulation. En revanche, le maclage n'apparaît pas dans la microstructure obtenue via l'élasticité non linéaire car le paysage énergétique du maclage n'a pas été pris en compte. Cette étape de validation nous a permis d'effectuer la même déformation mais cette fois sur un polycristal de TATB contenant 10 grains orientés aléatoirement. Les interactions élastiques entre grains génèrent des contraintes non-hydrostatiques très tôt et le flambage est déclenché pour la même compression que pour le monocristal (en DM et en élasticité non linéaire), affectant tous les grains. Les histogrammes contrainte-volume se dispersent très rapidement et évoluent alors principalement entre les les deux équations d'état présentées dans ce travail. Ces résultats représentent une étape importante dans la modélisation du comportement mécanique du monocristal et du polycristal de TATB. Les simulations en élasticité non linéaire reproduisent bien l'équation d'état en déformation sphérique obtenue en DM sur le monocristal et prédisent la microstructure de manière raisonnable. Les différences entre les équations d'état SD et HP devraient s'atténuer en présence d'un polycristal contenant de nombreux grains orientés aléatoirement. Avec un tel polycristal isotrope, la relation pression-volume devrait être la même sous pression hydrostatique et déformation sphérique. Cependant, le flambage de chaque grain du polycristal montre la complexité de son comportement et la nécessité, ainsi que le challenge, de considérer des méthodes d'homogénéisation compatibles avec cette pseudo transition de phase afin de reproduire une équation d'état représentative. During this thesis, we have investigated the mechanical behavior of the energetic material 1,3,5-triamino-2,4,6-trinitrobenzene using a multiscale approach. The only experimental data for TATB single crystal comes from X-ray diffraction and one nanoindentation testing performed recently. At the macroscopic scale, TATB-based plastic bonded explosives exhibit very defective TATB grains with porosity, fractures and large deformation bands that are believed to be the results of plastic activity. Since at this time no mechanical experiments can provide insight into the mechanical behavior at the microscopic scale, and regarding the necessity to build a constitutive law at mesoscopic scale, computational multiscale approach must be followed. It represents the core of this thesis. Thanks to the continuous improvement in high performance computing during the last decades, simulations at the scale where continuum mechanics can be used is at hand reach to improve our knowledge of the role of elastic and plastic anisotropy on the behavior of a milimeter size TATB polycrystal.

Fe 3

 3 afin d'étudier le comportement mécanique du TATB (1,3,5-triamino-2,4,6-trinitrobenzene). Les seules données expérimentales disponibles à l'échelle du monocristal proviennent d'expériences de diffraction aux rayons X et d'un essai de nanoindentation mené récemment. À l'échelle macroscopique, les explosifs comprimés au TATB présentent une structure parsemée de défauts tels que des pores, des fissures et de larges bandes de déformation que l'on suppose être la conséquence d'une activité plastique. Actuellement, très peu d'expériences sont à même de fournir des informations sur le comportement du TATB à l'échelle microscopique. De ce fait, la mise en place des lois de comportement aux échelles mésoscopiques et macroscopiques place l'approche multiéchelle au coeur de cette thèse. Grâce au progrès incessant qui a accompagné le calcul hautes performances ces dernières décennies, les simulations à l'échelle mésoscopique (µm) ont vu affiner leur résolution en temps et en espace afin de mieux tenir compte des détails des mécanismes de déformation à l'échelle atomique (nm). La dynamique moléculaire (DM) classique a été utilisée afin de simuler le monocristal de TATB pour diverses conditions de température et de pression. Ceci a permis le calcul de différentes propriétés physiques et mécaniques, permettant d'approfondir nos connaissances sur ce matériau très anisotrope pour lequel les résultats expérimentaux manquent cruellement. Dans ce travail, les simulations de DM (à l'échelle microscopique) ont été enrichies de diagnostiques spécifiques à la mécanique afin de permettre un dialogue avec l'échelle de la mécanique des milieux continus, étape inévitable si l'on veut construire des lois de comportement. Ces outils ont énormément aidé à comprendre le comportement du TATB à l'échelle microscopique mais se sont également avérés être très utiles lorsque des comparaisons entre les différentes techniques de simulation étaient nécessaires. Dans un premier temps, un potentiel inter-atomique provenant de la littérature a été utilisé afin de simuler le monocristal de TATB à l'équilibre via des simulations de DM avec des molécules rigides (chapitre 3). Les paramètres de maille ont été calculés via la méthode de Parrinello-Rahman pour des pressions allant jusqu'à 50 GPa et un bon accord avec l'expérience a été trouvé, même si le champ de forces prédit une distance inter-plans légèrement inférieure pour des conditions ambiantes. Partant des structures d'équilibre à différentes conditions de (T, P ), une nouvelle méthode pour calculer les constantes élastiques en température a été présentée et les résultats se sont avérés comparables avec ceux de la littérature. Une méthode similaire a également été présentée pour le calcul des γ-surfaces dans les plans basaux P 1 et P 2 du TATB en température et en pression. Dans les matériaux lamellaires, les dislocations sont connues pour nucléer et se déplacer facilement dans le plan de base (entre les plans atomiques/moléculaires). De ce fait, les structures de coeur des dislocations [100] et [010] dans les plans P 1 et P 2 du TATB ont été étudiées via le modèle de Peierls-Nabarro-Galerkin (PNG) dans le chapitre 4, représentant un premier pont multiéchelle. En raison des faibles énergies de faute d'empilement des plans basaux, les dislocations peuvent se diviser en plusieurs partielles, séparées par de larges bandes de faute d'empilement. La caractérisation de tels objets nécessiterait d'utiliser des grandes boîtes de simulations en DM. Pour contourner cette difficulté, les constantes élastiques ainsi que les γ-surfaces à 300 K à des pressions jusquà 10 GPa ont été utilisées dans le modèle PNG afin de calculer les structures de coeur des dislocations. Les dislocations (coin et vis) dans les plans P 1 et P 2 se séparent jusqu'à 4 partielles au coeur très complexe, étalées sur environ cent mailles de TATB. De plus, la pression n'a pas d'effet conséquent sur la structure de ces dislocations, en raison de la conservation de la proportion entre énergie élastique et inélastique sous pression. Les dislocations dans les plans basaux du TATB ne peuvent donc pas être considérées comme transverses au plan basal, les simulations de DM prédisent la nucléation de dislocations dans des plans transverses et quatre systèmes de glissement indépendants ont été identifiés. Ces dislocations sont nucléées via une rupture du réseau de liaisons hydrogènes des feuillets et possèdent un coeur très complexe comprenant une dilatance locale de l'ordre de 10 -15%. Les dislocations sont connues pour être le mécanisme principal contrôlant l'activité plastique dans les alliages métalliques et les minéraux. Dans l'optique de construire une loi de comportement en plasticité cristalline pour le TATB, ces premières informations seront à compléter par des données telles que la vitesse des dislocations, leurs modes d'interaction, leurs structures de coeur ou encore leurs contraintes d'activation. Ces données ne sont pas disponibles pour le TATB mais pourraient faire l'objet de travaux futurs. De plus, la majorité de ces propriétés ne pourrait pas être étudiée via la DM, en raison de la nécessité de considérer de très grandes boîtes de simulation afin de pouvoir prendre en compte les interactions élastiques à longues distance entre les dislocation par exemple. Cependant, l'effet de la pression sur les dislocations observées a révélé que la contrainte de nucléation était de plus en plus élevée. En effet, l'ajout d'une pression hydrostatique vient empêcher la dilatance locale du coeur, faisant ainsi monter la contrainte nécessaire à l'activation. En ce sens, il serait peut-être envisageable de ne considérer de tels mécanismes de plasticité que pour les sollicitations sous choc, par exemple. Lorsque la direction de compression se trouve dans les plans de base du TATB, une réponse très surprenante est observée. En effet, dès lors que la déformation longitudinale atteint ≈ 3%, une instabilité élastique de flambage est systématiquement déclenchée. Ce mécanisme, consistant en la succession de chevrons de déformation positive et négative, s'est avéré très dépendant de la direction de chargement. Les interfaces entre chevrons sont toujours perpendiculaires à la direction de compression et au plan basal. Plus la vitesse de déformation est faible, plus les chevrons sont larges. Une telle déformation peut être vue localement comme une composition entre un pur cisaillement dans le plan basal suivi d'une rotation. Afin de mieux comprendre la microstructure résultante, nous avons calculé une surface d'énergie en appliquant un cisaillement homogène dans le plan basal du monocristal de TATB. Ce calcul a mené à l'identification d'un chemin de déformation permettant de macler parfaitement la maille initiale de TATB. Nous avons cherché une éventuelle corrélation entre ce paysage énergétique de maclage et la microstructure observée lorsque le flambage est activé sous compression uniaxiale. Pour ce faire, nous avons extrait les distributions de la déformation des simulations de DM grâce à la mesure locale du tenseur gradient de la déformation. Les deux composantes non triviales de la déformation impliquée au cours du flambage et du cisaillement homogène ont été traitées puis reportées sur le paysage énergétique de maclage. Nous avons alors montré que le flambage dépend fortement de ce paysage énergétique. En effet, les distributions de la déformation se situent au voisinage des minimums du paysage énergétique. Lorsqu'une compression dans la direction [100] est effectuée, du flambage pur et équilibré est activé. Inversement, plus on comprime le TATB perpendiculairement à la direction [100], plus le maclage a tendance à s'activer. Ce nouveau mécanisme a été nommé mécanisme de Maclage-Flambage. Grâce à de précédents travaux sur le flambage de matériaux lamellaires et d'empilement de polymères, nous avons pu étudier la nucléation du flambage via la théorie de l'élasticité linéaire. Ainsi, si l'on décompose le flambage en une composition d'un cisaillement pur et d'une rotation d'ensemble, on peut lier les composantes de la déformation locale aux composantes de la déformation macroscopique appliquée au matériau, contenant une légère perturbation transverse à la direction de compression. Un seuil théorique de 4% a été trouvé, du même ordre de grandeur que celui observé en DM. Le fait que la nucléation de ce mécanisme soit donc purement élastique est intéressant et surprenant à la fois. En effet, cela signifie que dès lors qu'une compression uniaxiale supérieure à 3% est appliquée dans la direction des feuillets moléculaires du TATB, le mécanisme de Maclage-Flambage est activé et perturbe fortement la microstructure locale. Ceci implique également que ce mécanisme pourrait prévaloir sur les autres dans les premiers instants de la réponse mécanique du TATB. Une perspective intéressante serait de construire une loi de comportement qui prendrait en compte tous les mécanismes de déformation identifiés lors de ce travail de thèse. En revanche, puisque nous n'avons pas suffisamment d'informations concernant les dislocations non basales pour le moment, un modèle de plasticité cristalline n'est pas envisageable. Une solution serait d'adopter un modèle à seuil de type Norton-Bingham pour modéliser la plasticité dans les plans transverses du TATB ainsi qu'inclure une dépendance en pression du seuil afin de reproduire l'interaction de la dilatance avec la pression. D'autre part, le mécanisme d'activation de l'instabilité de flambage étant purement élastique, celui-ci est bien reproduit par l'élasticité non linéaire comme présenté dans le dernier chapitre et discuté ci-dessous. Enfin, une solution permettant de prendre en compte le chemin de transition menant au maclage serait d'utiliser la méthode des chemins de réaction en champ de phase utilisée par Denoual et al. [29] pour les transformations martensitiques dans un alliage modèle Ni ainsi que par Denoual [30] et Vattré [106] pour le comportement du fer sous haute pression. Dans le chapitre 6, nous nous sommes focalisés sur le calcul d'une équation d'état représentative pour le comportement polycristallin de matériaux très anisotropes avec une application sur le TATB. Deux équations d'état ont été considérées et ont mené à des résultats très différents. Lorsqu'une pression hydrostatique est appliquée sur le monocristal de TATB, chacun de ses paramètres de maille subit une évolution différente en accord avec la raideur correspondante et un champ de pression hy-drostatique est associé avec une déformation non sphérique. De manière opposée, lorsqu'une déformation sphérique est appliquée, chaque paramètre de maille subit la même déformation et une contrainte nonhydrostatique est induite. L'équation d'état du monocristal de TATB en pression hydrostatique (HP EOS) a été calculée via la DM et la méthode de Parrinello-Rahman tandis que l'équation d'état en déformation sphérique (SD EOS) a été calculée en utilisant la méthode des chemins de déformation, à partir de la même configuration de départ. Après une compression de 10%, une pseudo transition de phase correspondant au mécanisme de Maclage-Flambage est déclenchée lors du calcul de l'équation d'état SD EOS en déformation sphérique. En revanche, cette transition de phase du premier ordre n'apparaît pas en pression hydrostatique à cause de la grande stabilité du TATB lors de ce type de sollicitation. Ce résultat est très intéressant puisqu'il indique que sous pression, les matériaux fortement anisotropes ont une pseudo transition de phase induite uniquement par les contraintes non-hydrostatiques. Le comportement du monocristal de TATB en déformation sphérique a ensuite été simulé via un code Lagrangien de type "element free Galerkin" dans le contexte de la mécanique des milieux continus. Ce code utilise comme donnée d'entrée l'équation d'état SD EOS calculée en DM pour établir les constantes élastiques en pression (sous l'hypothèse de la conservation de l'anisotropie) pour une déformation sphérique jusqu'à V /V 0 = 0.85, suivie d'une relaxation à volume constant. La relation volume-pression obtenue est en accord qualitatif avec les résultats de l'équation d'état en déformation sphérique issus de la DM et l'instabilité de flambage est reproduite au même seuil de déformation. Après la nucléation de l'instabilité, la moyenne de l'évolution pression-volume dévie de la courbe SD EOS et évolue entre les courbes HP EOS et SD EOS considérées comme deux bornes. Cependant, le paysage énergétique de maclage n'ayant pas été pris en compte dans cette simulation, la microstructure obtenue a subi une déformation entièrement constituée de flambage, contrairement à la microstructure obtenue en DM qui consistait en un mélange de maclage et de flambage. La comparaison des niveaux de déformation a été rendue possible grâce au calcul local du tenseur gradient de la déformation en DM et des amplitudes similaires de la déformation ont été observées. Puisqu'un accord qualitatif entre les simulations de DM et les simulations en élasticité non linéaire du flambage a été observé, nous avons décidé de simuler un petit polycristal de TATB dans le même contexte. Dix grains avec des orientations aléatoires ont été générés via un pavage de Voronoï afin de construire un polycristal. Ensuite, une déformation sphérique macroscopique a été appliquée. En raison de l'orientation aléatoire de chaque grain, les interactions élastiques entre grains voisins provoquent des contraintes non-hydrostatiques et l'instabilité de flambage est déclenchée au même seuil que pour le monocristal. Cependant, la moyenne de l'évolution des histogrammes volume-pression dévie beaucoup plus tôt de la courbe SD EOS dans la simulation et continue d'évoluer entre les bornes HP EOS et SD EOS. Cette simulation a mené à une microstructure entièrement flambée pour chaque grain, montrant que la déformation sphérique macroscopique avait provoqué des contraintes non-hydrostatiques dans la totalité du polycristal. Ces différents résultats ont montré la complexité du comportement d'un polycristal possédant des grains très anisotropes comme le TATB. La nécessité de considérer des méthodes d'homogénéisation qui prendraient en compte cette pseudo transition de phase représente un réel challenge et devrait mener à des résultats très intéressants dans un futur proche.ism. The corresponding rotation matrix reads:

r 2 Figure A. 1 : 6 Figure B. 1 :B. 1

 21611 Figure A.1: Weight function w j (r) used for the computation of the local deformation gradient tensor with r normalized by the cutoff radius r cut used for the measure. In this work, r cut = 26 Å.

  Figure B.1(left). Multiple peaks of different intensity are clearly visible. However, due to the noisy property of the data, it is complicated to precisely identify local maxima positions in order to report them onto the twinning energy landscape presented in chapter 5 to assess possible correlations.

Figure B. 2 ( 6 Figure B. 2 :

 262 Figure B.2(left)) and labelled (see Figure B.2(rigth)) where the different colors represent the four different peaks found by this procedure. Finally, the different areas are labelled in order to apply the next operation that aims to compute the local maximum of each area.

Figure B. 3 (Figure B. 3 :Figure B. 4 :

 334 Figure B.3: Maxima positions and associated density. (legt) Gold stars represent the exact location of different maxima. (right) Gold stars have been replaced by circulaer shape of which the area is proportional to the number of molecules concerned by the associated deformation. The sum of all the molecules is more that 95%.

Keywords:

  TATB, multiscale modeling, plasticity, twinning-buckling, dislocations, anisotropy Modélisation Multiéchelle du Comportement des Matériaux Energétiques RESUME : La conception de lois de comportement en science des matériaux n'est pas nouvelle.Cependant, le progrès constant en calcul haute performance change la donne. En effet, ces lois visent désormais à tenir compte de la microstructure et de la physique sous-jacente, à l'échelle atomique, pour laquelle les techniques de simulation sont précises mais très coûteuses. L'approche multiéchelle semble parfaitement adaptée à ces problématiques et le dialogue entre échelles nécessaire. Dans cette thèse, le comportement mécanique du matériau énergétique TATB en température et en pression est étudié via des simulations de dynamique moléculaire afin de caractériser les mécanismes microscopiques responsable de son comportement irréversible. Le calcul local de variables mécaniques a été développé dans des simulations atomistiques, permettant le dialogue avec les méthodes continues. De plus, une méthode d'application de chemins de déformation a été couplée avec la dynamique moléculaire, menant à la caractérisation de la réponse mécanique très anisotrope du monocristal de TATB. Nucléation de dislocations au coeur complexe, chemin de transition pour le maclage et pseudo-transition de phase de type maclage-flambage sont trois comportements distincts associés à trois types de sollicitation dans différentes directions. Des simulations à l'échelle mésoscopique, alimentées par les données calculées à l'échelle microscopique, sont ensuite effectuées et visent à reproduire la pseudo-transition de phase sous compression triaxiale dans un code Lagrangien. La comparaison des résultats aux deux échelles est rendue possible par les outils de mécanique des milieux continus implémentés dans le code de dynamique moléculaire. Finalement, un polycristal de TATB est simulé en élasticité non linéaire et nous montrons l'importance de considérer une équation d'état compatible avec cette pseudo-transition de phase, qui semble avoir une forte influence sur le comportement du polycristal.Mots clés : TATB, modélisation multiéchelle, plasticité, maclage-flambage, dislocations, anisotropie

Table 2 .1: TATB lattice parameters at

 2

300 K, 0 GPa. a, b, c are in Å, α, β, γ in degrees.

Table 2 .2: Elastic Constants of TATB from literature at ambient condi- tions in

 2

	C 11	52.09	57.7
	C 22	50.18	58.0
	C 33	15.75	17.0
	C 44	1.20	1.0
	C 55	0.73	0.6
	C 66	17.58	20.3
	C 12	14.24	16.2
	C 13	2.05	3.2
	C 14	-0.14	0.1
	C 15	-0.61	-0.9
	C 16	0.30	0.0
	C 23	3.55	5.7
	C 24	0.53	0.6
	C 25	-0.15	-0.5
	C 26	0.91	2.0
	C 34	0.05	-0.1
	C 35	-0.28	-0.3
	C 36	-0.14	-1.0
	C 45	0.07	0.01
	C 46	-0.16	-0.5
	C 56	0.20	0.1

Bedrov et al.

including covalent bond stretching. They also reported values at 0 K through MD simulations by least-squares fitting of strain-energy-density vs. strain curves. Both predictions of elastic constants at ambient conditions are indicated in Table

2

.2 using the conventional Voigt notation. It can be noted that the results obtained with the force field including covalent bond stretching predict a softer elastic response. Furthermore, the elastic tensor has many components near zero and could be assimilated to the stiffness matrix of orthotropic materials which have 9 independent elastic constants and a very soft shear behavior in the basal plane.

The understanding of structural defects and irreversible mecha-Elastic Constants MD-MC [1] MD-MC [2] GPa. [1] [73], [2] [11]

2.4.1 Structure et propriétés du TATB monocristallin 2.4.1.1 Structure cristalline

	supplémentaires ont été nombreux [54, 98, 39, 23, 88]. En 2007, Stevens
	Le présent chapitre est consacré à un bref examen des travaux exis-et al. [100] ont déterminé une équation d'état du TATB jusqu'à 13 GPa
	tants sur le comportement du TATB. Nous rappelons en partie 2.4.1 les via une expérience de diffraction aux rayons X en enclume de diamant.
	différentes études qui ont été réalisées sur le monocristal, l'accent étant Plus récemment, Plisson et al. [91] ont effectué une expérience similaire
	mis sur les quelques résultats expérimentaux disponibles ainsi que sur la jusqu'à 66 GPa. Différents résultats expérimentaux sont reportés dans le
	présentation des modèles atomistiques développés pour le TATB ayant Tableau 2.1.
	mené à l'obtention de propriétés thermo-mécaniques via des simulations
	de DM. Nous présentons ensuite de manière succinte des résultats à la 2.4
	fois numériques et expérimentaux concernant le comportement macro-
	scopique de matériaux composés à base de TATB en partie 2.4.2.
	Le TATB est un matériau énergétique utilisé dans de nombreuses
	formulations d'explosifs. Il possède une maille triclinique de groupe de
	symétrie P1 [19] contenant deux molécules de formule C 6 H 6 N 6 O 6 (voir
	Figure 1.1). Ces molécules consistent en un noyau benzénique ayant
	une alternance de groupes (NO 2 ) et (NH 2 ), avec de fortes liaisons intra-
	moléculaires. Ce matériau est référencé comme un "matériau 2D de van
	der Waals" à cause du fort contraste entre les liaisons hydrogène inter-
	moléculaires au sein d'un même feuillet et les interactions entre molécules
	appartenant à différents plans principalement gouvernées par des inter-
	actions faibles de van der Waals. Ces différents types d'interactions, ainsi
	que l'empilement graphitique des molécules, confèrent au TATB des pro-
	priétés thermiques et mécaniques très anisotropes [11, 58, 59].
	Les premiers travaux sur le cristal de TATB ont été menés par Cady
	et Larson [19] en 1965 via une expérience de diffraction aux rayons X sur
	un échantillon de 200 µm. Ils ont été les premiers à identifier la structure
	du TATB ainsi que le réseau de fortes liaisons hydrogène lui conférant
	une grande stabilité sous pression. Depuis, les travaux visant à mieux
	caractériser la structure du TATB et à trouver de potentielles phases

.1.2 Potentiels atomistiques pour le TATB

	Plusieurs champs de forces ont été proposés pour modéliser le
	TATB via des simulations de DM [44, 94, 11, 58]. Ces champs de forces
	peuvent être réactifs et traiter les molécules comme rigides. Cette dernière
	considération exclue toute flexibilité comme avec le potentiel de Rai [94]
	par exemple. Gee et al. [44] ont quant à eux développé un potentiel pour
	molécules flexibles, autorisant tout type de distortion. Finalement, Bedrov
	et al. [11] ont construit un champ de forces (deux versions, une polarisée,
	l'autre non) qui tient compte des valeurs à l'équilibre des liaisons intra-
	moléculaires. Puisque la version non polarisée de ce dernier a été utilisée
	lors de cette thèse, il sera détaillé au prochain chapitre.
	Les champs de forces énumérés ci-dessus ont tous été construits
	pour reproduidre la structure du TATB et certains sont paramétrés via des
	simulations de chimie quantique et autres données expérimentales. La
	principale différence entre ces modèles provient de la forme utilisée pour
	les interactions inter-atomiques : Lennard-Jones pour Gee et al. et Rai et
	al., Buckingham (exp -6) pour Bedrov et al.. En raison du manque de résultats expérimentaux concernant les
	propriétés thermo-mécaniques du TATB monocristallin, les simulations
	de DM ont été largement utilisées ces dernières années. Un brève revue
	de ces études est proposée dans ce qui suit.

2.4.1.3 Propriétés mécaniques

	jusque là. Bedrov et al. [11] ont calculé le tenseur d'élasticité en utilisant calculer la dureté du plan (001) avec une valeur de 1, 02 GPa à l'aide d'un
	la méthode des fluctuations de Parrinello-Rahman [86] via une approche indenteur de type Berkovich. Leurs résultats sont en très bon accord avec
	de Monte Carlo (DM-MC) [16, 97]. Mathew et al. [73] ont également cal-l'expérience menée récemment par Taw et al. [102].
	culé ce tenseur avec une approche similaire mais en utilisant une version
	modifiée du potentiel [58] de Bedrov et al. Les valeurs des constantes 2.4.1.4 Propriétés thermiques
	élastiques obtenues au cours de ces deux études sont données dans le Plusieurs propriétés thermiques ont été étudiées ces dernières an-
	Tableau 2.2. nées. Kroonblawd et Sewell ont calculé la conductivité thermique du TATB
	La compréhension de la déformation irréversible du TATB reste li-selon trois directions [58], celles dans les feuillets de molécules étant 1.7
	mitée. De manière générale, les dislocations interviennent au premier fois plus grandes que celle dans la direction transverse. Ce travail fut
	plan lorsqu'il s'agit de la plasticité des matériaux. Pour les matériaux la-étendu à l'étude des effets de taille et de la présence de lacunes sur la
	mellaires, Delavignette [25] a montré que la déformation plastique est conductivité thermique [59].
	essentiellement gouvernée par le glissement de dislocations dans des Plus récemment, Kroonblawd et Sewell ont étudié la relaxation ani-
	plans parallèles aux feuillets (appelés plans de base). Bien que quelques sotrope de points chauds dans le TATB pour des échelles allant jusqu'à
	exceptions existent, ces dislocations sont souvent sujettes à une dissocia-tion en deux ou plusieurs partielles [25, 4]. Cette dissociation est connue ≈ 10 nm via des simulations de DM [60]. Ils ont considéré deux direc-tions de conductivité, en accord avec l'anisotropie observée dans leurs
	pour être contrôlée par les énergies de faute d'empilement (γ SF ) le long de travaux précédents et ont obtenu un très bon accord avec les prédictions
	la direction de glissement. Afin d'étudier ce phénomène, Mathew et al. [73] analytiques.
	ont calculé les surfaces d'énergie de faute d'empilement (γ-surfaces) dans Kroonblawd et al. ont étudié la cinétique et la dynamique d'échange
	les plans de base (001). Deux plans de base doivent être considérés du d'énergie entre les degrés de libertés internes et externes via des simula-
	fait de la centrosymétrie du matériau. De faibles valeurs pour les énergies tions de DM afin d'adapter un modèle de coarse-graining en dynamique
	de fautes d'empilement ont été calculées, suggérant un glissement très des particules dissipatives (DPDE) [56]. La courbe de fusion ainsi que des
	facile des dislocations ainsi qu'une contrainte très faible nécessaire à leur propriétés du TATB liquide ont été récemment calculées en température
	mobilité. et en pression par Mathew et al. [70] qui ont utilisé des simulations de DM
	Mathew et al. [74] ont présenté des simulations de nanoindenta-avec un potentiel non réactif et des molécules flexibles.
	tion via des simulations de DM à 77 K sur les plans (100) (010) et (001) en
	utilisant un indenteur rigide et sphérique. Ils ont observé une réponse
	hertzienne sur le plan (001) avec une déformation inélastique via de la
	délamination et de l'ondulation des feuillets. Une réponse non-hertzienne
	fut observée pour les deux autres plans. Des modules d'indentation de
	27, 7, 26, 0 and 10, 7 GPa pour les plans (100), (010) et (001) ont été prédits,
	Aucune mesure des constantes élastiques du TATB n'a pu être réalisée traduisant la forte anisotropie du TATB. Les auteurs ont également pu

2.4.2 Propriétés d'un composé explosif au TATB 2.4.2.1 Microstructure

	granulée en billes millimétriques. Après une étape de séchage, un pres-
	sage isostatique sous vide et en bain d'huile est effectué à température
	modérée et haute pression. Une composition au TATB à faible porosité
	(< 2, 5%) est alors obtenue, avec de bonnes propriétés mécaniques iso-
	tropes et une aptitude à l'usinage [66, 3]. La microstructure finale, avec
	des grains de l'ordre de 50 µm, est présentée en Figure 2.1. Une porosité
	résiduelle peut être observée sous la forme de petits points noirs tan-
	dis que d'autres défauts "mécaniques" sont très facilement identifiables :
	des fissures intra-granulaires sont présentes (ellipses rouges) ainsi que
	des larges bandes de déformation (cercle vert), non présentes avant le
	pressage isostatique, qui sont certainement le fruit d'une déformation
	plastique.
	2
	À l'échelle macroscopique, des aggrégats polycristallins de TATB
	sont associés avec un liant polymérique, dont la fraction est inférieure
	à 5% en volume. Le processus de fabrication comporte deux étapes : la
	poudre de TATB est enrobée d'un liant en phase liquide tout en étant

.4.2.2 Comportement mécanique

  GPa et 21, 2 GPa, équivalent à un coefficient de Poisson de 0, 335.

	de 7, 1 Ceci permet l'estimation du module de compressibilité (7, 2 GPa) et de
	cisaillement (2, 7 GPa). De plus, un coefficient de dilatation thermique
	volumétrique de 1, 62 × 10 -4 K -1 a été trouvé, environ deux fois plus petit que celui du monocristal de TATB présenté en partie 2.4.1.4.
	De nombreux travaux en cours s'attachent à simuler le comporte-
	ment macroscopique des matériaux composés au TATB. Ambos et al. [3]
	ont proposé une modélisation numérique de l'expansion thermique d'un
	explosif au TATB, tandis que Gasnier et al. [43] ont calculé ses propriétés
	thermo-élastiques par une modélisation en champ complet de type FFT.
	Récemment Gasnier [42] a commencé a étudier l'impact de la présence
	de fissures sur les propriétés macroscopiques. Cependant, ces modèles
	ne prennent pas en considération les mécanismes de déformation du mo-
	nocristal de TATB puisqu'aucune donnée (expérimentale et numérique)
	Le comportement quasi-statique de cette composition de TATB a été n'est actuellement disponible. Le but principal de cette thèse est donc
	largement caractérisée via un grand spectre de mesures expérimentales : d'étudier les mécanismes de déformation du TATB monocristallin afin
	compression simple, monotone, cyclique, compression tri-axiale à basse de mieux comprendre ce matériau très complexe et d'alimenter les tech-
	et haute pression etc. [3]. En traction simple, le comportement fragile [15, niques de modélisation aux échelles supérieures.
	103] provoque une rupture dans le domaine élastique tandis qu'en com-
	pression simple l'endommagement est progressif jusqu'à la rupture	à
	environ 2% de déformation longitudinale. Une évolution en pression
	de la déformation à rupture est également observée. Une déformation
	longitudinale de 1, 5% pour un test de compression simple mène à une
	contrainte d'écoulement de 30 MPa tandis que pour un pressage confiné
	de 200 MPa, la contrainte seuil approche les 85 MPa pour une déformation
	aux alentours de 30%. Ce comportement est typique d'une transition
	fragile-ductile observée dans les matériaux quasi-fragiles sous pression.
	2.4.2.3 Comportement thermoélastique
	À l'échelle du polycristal, des modules élastiques longitudinaux
	et transverses ont été mesurés expérimentalement [3] avec des valeurs

2.4.3 Conclusion

	comportement aux échelles supérieures. 3
	CHAPTER 3
	COMPUTATIONAL METHODS AND FIRST VALIDATION
	RESULTS
	The main computational methods used in this work are presented. Sec-
	tion 3.1 concerns the classical molecular dynamics (MD) method, in-
	cluding a detailed description of the force field used to represent TATB
	microscopic behavior. Furthermore, results on the elastic constants in
	temperature and pressure as well as the equation of state are presented
	as validation studies. Then, section 3.2 explains the methodology used
	to prescribe deformation paths in MD as well as the computation of a
	local deformation gradient tensor in a Lagrangian formulation that links
	microscopic observations during MD simulations to what can be obtained
	in continuum mechanics. Finally, section 3.3 is dedicated to the intro-
	duction of the Peierls-Nabarro-Galerkin (PNG) model, used in this work
	to compute the core structure of dislocations in the basal (001) plane of
	TATB single crystal.
	Parts of this chapter previously appeared in:
	La revue de la littérature présentée ci-dessus a montré la difficulté [62], Dislocation Core Structure at Finite Temperature Inferred by Molecular Dy-
	que représente l'obtention de données expérimentales sur le comporte-namics Simulations for 1,3,5-Triamino-2,4,6-Trinitrobenzene Single Crystal, J. Phys.
	ment mécanique du TATB. En effet, bien que la structure cristalline du Chem. C, vol. 121, no. 23, pp. 235901, 2017
	monocristal de TATB soit connue via la diffraction aux rayons X, une seule [63], Irreversible Deformation Mechanisms for 1,3,5-Triamino-2,4,6-Trinitrobenzene
	expérience mécanique de nanoindentation à l'échelle du grain a pu être Single Crystal Through Molecular Dynamics Simulations, J. Phys. Chem. C, In edition,
	réalisée récemment par Taw et al. [102]. Quelques études en simulation 2018
	numériques ont été effectuées par Mathew et al. qui ont étudié les pro-
	priétés du glissement dans le plan basal du TATB [73] et sa réponse sous
	nanoindentation [74]. Le but principal de cette thèse est donc d'explorer
	les éventuels mécanismes de déformation du TATB monocristallin afin de
	garantir des fondations solides pour la construction d'une future loi de
	39

.1 Classical molecular dynamics 3.1.1 First principles

  .1 concerning bending

		Atomic	Atomic
		polarizability	charge	A	B	C
	Atom type	α( Å3 )	q(e)	(kcal mol -1 )	( Å-1 )	(kcal Å3 mol -1 )
	C(-NO 2 )	0.70	0.2839	107 023.9	3.6405	554.01
	C(-NH 2 )	0.70	-0.2025	107 023.9	3.6405	554.01
	N(-O 2 )	1.00	0.5665	13 783.9	3.3333	423.36
	N(-H 2 )	1.00	-0.5775	39 091.8	3.3158	833.48
	O	1.00	-0.4252	15 923.1	3.6446	239.07
	H	0.21	0.3900	7 584.2	5.2846	8.23
		Bond length	Bend type

Table 3 .1: TATB force field parameters from Bedrov et al. [11]

 3

  This work Uncertainty (σ) Deviation (%)

	a( Å) b( Å) c( Å) α( • ) β( • ) γ( • )	9.01 9.03 6.81 108.58 91.82 119.97	9.01 9.04 6.57 103.34 88.34 120.70	1.2 × 10 -3 9.1 × 10 -4 5.9 × 10 -4 7.8 × 10 -3 1.5 × 10 -2 6.8 × 10 -3	-0.01 0.14 -3.55 -4.83 -3.79 0.61

Table 3 .2: Unit cell parameters from MD Parrinello-Rahman simula- tions compared to experimental values at ambient condi- tions

 3

. are compared to the two experiments driven by Stevens et al. [100] and Olinger and Cady [80] respectively, and with a density functionnal theory study of Pineau et al.

Table 3 .3: Variation with T of the TATB elastic coefficients at

 3 .4 are the stiffnesses in reference state (C ijkl ).

		0 K	200 K 250 K 300 K 350 K 400 K
	C	56.10 51.72 51.45 49.88 48.96 49.16
	C	57.38 52.83 51.53 50.07 49.94 48.27
	C	26.77 24.94 24.83 24.66 24.18 23.66
	C	3.58	2.64	2.67	2.32	2.22	2.27
	C	3.51	2.04	2.08	1.92	1.86	1.83
	C	30.42 27.14 26.40 25.22 24.92 24.22
	C	11.69 11.83 11.83 11.41 11.67 11.79
	C	2.52	2.73	3.04	2.39	2.34	2.74
	C	-1.01 -0.11 -0.06 -0.24 -0.23 -0.14
	C	-2.10 -1.50 -1.15 -1.52 -1.42 -0.79
	C	1.21	2.52	2.72	2.47	2.32	2.86
	C	5.80	4.49	4.48	4.30	4.09	3.99
	C	-1.97 -0.52 -0.46 -0.75 -0.64 -0.45
	C	1.15 -0.09 -0.07 0.29	0.33	0.03
	C	4.31	3.14	3.05	3.75	3.75	3.75
	C	-1.15 -1.27 -1.09 -1.37 -1.17 -0.78
	C	-0.24 -0.85 -0.88 -0.80 -0.77 -0.84
	C	-0.66 -0.62 -0.53 -0.79 -0.69 -0.38
	C	0.06 -0.20 -0.20 -0.07 -0.07 -0.08
	C	0.92 -0.09 -0.08 0.29	0.32	0.02
	C	-0.94 0.12	0.16 -0.14 -0.10 0.07

0 GPa. Values are in GPa.

2

  GPa 4 GPa 6 GPa 8 GPa 10 GPa

	C	126.17 184.13 234.56 278.77 325.38
	C	108.00 161.96 213.81 265.28 308.08
	C	53.92	72.55	91.16 110.00 128.61
	C	7.65	11.97	17.21	21.90	25.74
	C	7.52	12.75	16.31	20.42	24.31
	C	54.18	79.73 103.30 126.90 147.67
	C	27.22	42.42	57.45	72.04	84.67
	C	6.95	12.38	17.17	21.68	25.93
	C	-1.58	-2.04	-1.82	-1.87	-2.12
	C	-1.08	-0.19	0.71	0.69	-0.32
	C	2.85	2.11	1.55	0.67	-0.58
	C	10.81	15.67	22.75	29.45	34.07
	C	-5.18	-7.59	-7.37	-7.72	-9.34
	C	0.73	1.19	1.80	1.98	2.33
	C	10.58	14.33	17.04	22.08	27.26
	C	-2.11	-2.96	-1.68	-0.60	0.06
	C	-0.34	1.19	2.21	1.41	1.07
	C	-2.12	-3.92	-5.05	-6.12	-6.57
	C	0.32	-0.36	-0.44	-0.44	0.36
	C	0.74	1.19	1.79	1.94	2.29
	C	-0.79	-1.76	-2.69	-3.89	-4.99

Table 3 .

 3

4: Variation with P of the TATB elastic coefficients at 300 K. Values are in GPa.

  Afin de calculer les forces entre particules, les potentiels sont en géneral tronqués pour une distance égale à r cut appelée le rayon de coupure. Des corrections sont ensuite appliquées afin d'assurer la continuité des forces et de l'énergie à cette distance. Cette troncature réduit ainsi le nombre de voisins avec lesquels une particule va devoir interagir. On considère une liste de voisins appelée liste de Verlet[START_REF] Verlet | Computer experiments on classical fluids. ii. equilibrium correlation functions[END_REF] dans une sphère environnante de rayon r v > r cut . Cette liste reste inchangée tant que le déplacement d'une particule ne dépasse pas r v -r cut , ce qui représente une économie de calcul conséquente.

	suit :	1. v(t + 2. r(t + ∆t) = r(t) + v(t + 1 2 ∆t) = v(t) + 1 a(t)∆t 2 1 3. v(t + ∆t) = v(t + 1 2 ) + 2 a(t + ∆t)∆t 1 2 ∆t)∆t	(3.29)
	Ensembles thermodynamiques
	Ensemble microcanonique
	L'ensemble microcanonique permet d'étudier un système isolé compre-

.28) o ù V est la vitesse du CDM, m est la masse, I est le tenseur d'inertie autour d'axes passant par le CDM, ω la vitesse angulaire et U l'énergie potentielle.

Intégration numérique des équations du mouvement

Afin de calculer les positions des particules à l'instant t + ∆t, l'algorithme de Verlet-Vitesse

[START_REF] Verlet | Computer experiments on classical fluids. ii. equilibrium correlation functions[END_REF]

a été utilisé. Il consiste à évaluer dans un premier temps les vitesses à t + ∆t 2 avant de calculer r et v à t + ∆t comme L'étude des trajectoires atomiques en dynamique moléculaire permet l'obtention de certaines propriétés "macroscopiques". Deux ensembles thermodynamiques ont été utilisés dans ce travail et sont dénotés NVE et NVT [105, 31]. Un autre ensemble, celui de Parrinello-Rahman, qui comprend une modification du Lagrangien a été utilisé et est décrit en annexe A.2. nant N particules dans un volume V d'énergie totale E. Les simulations se font donc à N , V et E constants puisqu'aucun échange n'est possible avec l'extérieur. La fonction de partition Q NVE permet de décrire la distribution et la probabilité des états du système :

3.5.2.2 Mesure locale du tenseur gradient de la transformation

  sur les zones du matériau reconstruites après le passage d'un défaut, ou encore sur des zones qui ont uniquement subit des grandes rotations. Pour ces raisons, le calcul local du tenseur gradient de la transformation F est présenté dans ce manuscrit. Cette mesure a déjà été introduite par plusieurs auteurs (Gullett et al.[START_REF] Gullett | A deformation gradient tensor and strain tensors for atomistic simulations[END_REF], Tucker et al.[START_REF] Tucker | Continuum metrics for deformation and microrotation from atomistic simulations: Application to grain boundaries[END_REF], Denoual et al.[START_REF] Denoual | Phase-field reaction-pathway kinetics of martensitic transformations in a model fe 3 Ni alloy[END_REF]) et plus récemment par Kober et al.[START_REF] Kober | Characterizing atomistic geometries and potential functions using strain functionals[END_REF].

	En mécanique des milieux continus, un point du matériau au temps
	t dans la configuration déformée est exprimé en fonction de sa position
	dans l'état de référence via le tenseur gradient de la transformation par

La déformation appliquée au cours des simulations de DM peut mener à la nucléation de défauts locaux au sein du cristal. Deux diagnostics peuvent être appliqués pour identifier ces défauts. Premièrement, lorsque des dislocations sont impliquées, le cristal est reconstruit mais fortement perturbé au voisinage du coeur de la dislocation. Une détection locale de la modification structurelle permet d'identifier la ligne de dislocation, mais pas le plan. Pour des défauts comme le maclage, seuls les plans de macle sont identifiés. Un tel diagnostic tient compte uniquement de l'état actuel du cristal sans prendre en considération l'histoire de la déformation. En revanche, lors d'une mesure Lagrangienne de la déformation, qui consiste à mesurer l'écart entre la configuration actuelle et la configuration initiale, toute l'histoire de la déformation est capturée. Le tenseur gradient de la transformation local F , décomposé en une rotation et une pure élongation (décomposition polaire F = R • Û ) fournit des indications essentielles x(t) = F (t) • X, o ù x et X sont les positions courantes et initiales. On approche le tenseur gradient de la transformation (appliqué lors de la simulation) par un tenseur local F i , défini pour chaque molécule par une méthode des moindres carrés qui se base sur les déplacements relatifs des molécules contenues dans le voisinage. La définition du formalisme utilisé est donnée dans l'annexe A.3.

Table 5 .2: Size effect on the nucleation stress

 5

	75 90	1.618 2.448 2.269	1.346 2.097 1.834	≃ 20.2 ≃ 16.7 ≃ 23.7

σ v (θ, φ) for θ = 90 •

Table 5 .3: Strain rate effect on the nucleation stress

 5

	15 30 60 75 90	612 1.409 1.249 1.699 1.781 1.833	1.407 1.215 1.082 1.458 1.537 1.593	≃ 12.7 ≃ 13.8 ≃ 13.4 ≃ 14.2 ≃ 13.7 ≃ 13.1

σ v (θ, φ) for fixed θ = 0 •

  . Le module de raideur K 0 et sa dérivée K ′ 0 sont identifiés afin d'approcher les données expérimentales ou numériques. De telles équations d'état dépendent de la variation de volume, du module de raideur et de sa dérivée mais ne précisent pas si la densification provient de l'application d'une contrainte hydrostatique ou non-hydrostatique. Lorsqu'un monocristal de TATB est soumis à une contrainte hydrostatique, son anisotropie élastique est conservée et chaque paramètre de maille suit une évolution différente, cohérente avec la valeur de la constante élastique associée. De ce fait, un volume sphérique initial se verra transformé en une ellipsode soumise à un champ de contrainte parfaitement hydrostatique (voir Figure 6.1). Cet état de contrainte mène à l'obtention du module de raideur dépendant d'une combinaison des souplesses S ij , obtenu via le tenseur des souplesses élastiques (inverse du tenseur de raideurs élastiques) et l'équation d'état associée est nommée "HP EOS" pour "hydrostatic pressure equation of state".

	transmetteur se solidifie au delà d'une certaine pression, des contraintes
	non-hydrostatiques sont générées. En présence de poudre, les conditions
	hydrostatiques sont assurées uniquement si celle-ci est bien dispersée
	dans le liquide, contrairement au cas o ù des contacts entre grains auraient
	lieu.
	Nous présentons par la suite des résultats de dynamique moléculaire
	concernant les deux équations d'état discutées précédemment ainsi qu'une
	première modélisation en mécanique des milieux continus à l'aide d'un
	code Lagrangien.
	Lorsque le même monocristal est soumis à une déformation sphéri-
	que imposée, o ù chaque direction suit la même évolution, la pression
	n'est plus contrôlée. De ce fait, le volume sphérique initial subit une
	simple homothétie et l'élasticité anisotrope du TATB mène à un champ
	de contrainte non-hydrostatique (voir Figure 6.1). Ici, l'hypothèse sur les
	déformations mène à l'obtention du module de raideur dépendant d'une
	combinaison des C ij et l'équation d'état associée est nommée "SD EOS"
	pour "spherical deformation equation of state". Elle est différente de HP
	EOS en général.
	Les conditions en pression hydrostatique sont souvent recherchées
	expérimentalement. Cependant, ces conditions sont souvent difficiles	à
	garantir à cause du matériau mis en jeu ou du milieu transmetteur de
	pression utilisé, et peuvent mener à des résultats discutables pour des
	matériaux très anisotropes comme le TATB. Si un monocristal parfait
	entouré d'un milieu transmetteur liquide est dans la cellule à enclumes de
	diamant, une pression hydrostatique est garantie. A contrario, si le liquide

6.7.2 Calcul d'équations d'état en dynamique moléculaire

La méthode de Parrinello-Rahman permet d'appliquer une consigne de pression hydrostatique sur le système simulé. Cette technique (présentée en annexe A.2) consiste à minimiser l'énergie relativement aux paramètres de maille afin d'assurer une contrainte hydrostatique. On a donc :

1.2.2 Cristaux moléculaires énergétiquesDans le vaste domaine des matériaux, les cristaux moléculaires occupent une place bien spéciale. De manière analogue aux céramiques ou encore aux métaux, ils peuvent être caractérisés par leur structure cristalline et par la position des molécules au sein de leur maille primitive.En effet, ces matériaux consistent en général en un arrangement régulier de molécules complexes, définissant ainsi un empilement périodique de basse symétrie. Les cristaux moléculaires, utilisés dans de nombreux domaines, par exemple par les entreprises pharmaceutiques[START_REF] Wildfong | Towards an understanding of the structurally based potential for mechanically activated disordering of small molecule organic crystals[END_REF][START_REF] Bates | Assessment of defects and amorphous structure produced in raffinose pentahydrate upon dehydration[END_REF][START_REF] Feng | Toward understanding the evolution of griseofulvin crystal structure to a mesophase after cryogenic milling[END_REF] et les compositions d'explosifs, ont été très largement étudiés au cours des dernières décennies. Par exemple, Lei et al.[START_REF] Lei | Mesoscale modeling of dislocations in molecular crystals[END_REF] ont étudié l'évolution des dislocations dans les cristaux moléculaires du sucrose et du paracétamol, tandis que les matériaux énergétiques (ME) tels que le cyclotriméthylènetrinitriamine (RDX), le tétranitrate de pentaérythritol (PETN), le cyclotétramé-

Cf. p161

Cf. p159

des lignes discrètes mais plutôt comme des larges bandes de faute d'empilement pour les conditions étudiées. La capacité du TATB à glisser très facilement dans ces directions soulève la question de l'association d'un paramètre de viscosité à ce glissement afin de le prendre tout de même en compte à l'échelle mésoscopique. Afin d'apporter un point de vue mécanique aux simulations de DM, nous avons implémenté une méthode qui permet d'imposer un chemin de déformation quelconque à la boîte de simulation de DM. Ainsi, on peut appliquer des cisaillements directionnels ou des compressions uniaxiales et triaxiales, de manière similaire à ce qui est fait dans des codes de mécanique des milieux continus Lagrangiens. Ces grandes déformations provoquent souvent la nucléation de défauts tels que des dislocations ou du maclage, entraînant des fortes perturbations locales et des rotations d'ensemble. Nous voulions donc être capables d'identifier et de caractériser ces phénomènes lors des simulations de DM. Afin de répondre à ce besoin, nous avons implémenté dans le code de DM une méthode permettant de calculer le tenseur gradient de la déformation en chaque centre de masse des atomes/molécules. Nous avons alors facilement accès à l'histoire complète du matériau via une décomposition polaire. Nous pouvons alors identifier les plans ayant subi le passage de dislocations ou encore les zones de fortes rotations, signatures typiques du maclage. La méthode des chemins de déformations en DM a été utilisée dans le chapitre 5 afin d'appliquer des cisaillements directionnels à un monocristal de TATB . Ce cisaillement est construit à partir d'un vecteur de compression et d'un vecteur de dilatation dans une direction transverse. Une surface de contrainte seuil de nucléation en cisaillement a été calculée et a permis l'identification de trois comportements distincts. Lorsque le cisaillement est appliqué dans le plan basal du TATB, un glissement homogène est activé, rejoignant ainsi les conclusions sur les structures de dislocation dans ces plans là. En revanche, lorsque le cisaillement est appliqué dans des directions

être considérées comme des lignes discrètes mais plutôt comme de larges bandes de faute d'empilement, pour les pressions et les températures étudiées (0 -10 GPa et 0 -300 K).

CHAPTER 5 MECHANISMS OF THE IRREVERSIBLE DEFORMATION

Introduction

The plastic behavior of the insensitive energetic molecular crystal TATB under moderate to large strains is investigated through MD simulations in this chapter. The method presented in section 3.2.1, designed to prescribe any deformation path, is used to apply various deformations to a TATB single crystal. In chapter 4, the plastic activity in the (001) TATB basal planes has been investigated by characterizing the dislocation core structures for pressure up to 10 GPa. Due to the splitting of dislocations into up to 4 partials, spreaded over a length approximately equals to a hundred lattices, dislocations cannot be considered as discrete lines but rather as large stacking fault ribbons.

In order to build a constitutive law this information has to be completed by a more comprehensive loading through a simulation set-up dedicated to the analysis of the onset and propagation of irreversible deformation mechanisms in TATB single crystal. Stress-strain curves are obtained for common deformations such as pure shear and uniaxial compressions, with loading directions spanning all possible orientations. Furthermore, using the Lagrangian local estimate of the deformation gradient tensor presented in section 3.2.2 and explained in detail in appendix A. For a deformation consisting in a compression along the y axis (ε yy ≤ 0 and εzz ≥ 0), a decrease is induced in the squared local deformations, leading to the generation of an instability. With the expressions of Equations 5.6 and 5.10, the averaged elastic energy reads:

where = 0 and thus the potential limit at which buckling occurs. If pure compression is applied (i.e. with no stretch along the normal to the layers), buckling is shown to happen at εyy = -0.043 (represented as a star). The effect of the transverse stretch εzz is clearly visible in Figure 5.17: compression along the z axis will prevent the apparition of buckling, i.e., the value of εyy needed to reach the instability becomes larger. Conversely, a tensile deformation along the z axis will promote the apparition of buckling while lowering the threshold value of εyy . For example, with εzz = 0.05, the compression needed to buckle falls from εyy = -0.043 to εyy = -0.022 which indicates the strong influence of the transverse deformation for buckling instability.

Compared to the threshold at which buckling occurs in MD simulations (around 3 %), the values of 4.3 % is slightly greater. However, this estimation is based only on the components of the TATB elastic tensor and particularly on the value of C 55 which has been measured through Olinger and Cady (grey disks) deviates from the other results, certainly a consequence of the two assumptions made by the authors on the crystal structure. However, we remark that their data are following the HP EOS which is consistent with the hydrostatic pressure conditions targeted by their experiment.

More is to say on the data of Stevens et al.: the blue disks follow well the HP EOS even if above 10 GPa a perturbation appears in the evolution, which may be a signature of the freezing of the ethanol/methanol solution above this pressure. They used three semi-empirical EOS (Murnaghan [START_REF] Murnaghan | The compressibility of media under extreme pressures[END_REF], Vinet [START_REF] Vinet | A universal equation of state for solids[END_REF] and Birch-Murnaghan [START_REF] Murnaghan | Finite deformations of an elastic solid[END_REF][START_REF] Murnaghan | The compressibility of media under extreme pressures[END_REF][START_REF] Birch | Finite elastic strain of cubic crystals[END_REF]) to fit the data. With the Vinet EOS, they reported values of 17.5 GPa and 7.6 for K 0 and K ′ 0 over the entire P -V range.

The data of Plisson et al. is represented as green disks in Figure 6.5 with a good accordance with the HP EOS curve in dark red. The authors used a Vinet EOS to fit over the whole experimental data leading to K 0 = 24.9 GPa, K ′ 0 = 7.4 and a smaller initial volume than the experimental data. Since a single Vinet EOS failed at properly fitting over the entire range of pressure, they splitted their data in two regimes, with a separation at P = 10 GPa. This resulted in K 0 = 7.2 GPa and K ′ 0 = 13.9 under 10 GPa and K 0 = 36.8 GPa and K ′ 0 = 6.6 for the remaining part. The authors attributed this two regimes evolution to the anisotropic compression ratio that present a crossover between 5 and 20 GPa. However, the presence of multiple single crystals in the DAC may have led to interaction between different grains and thus apparition of non-hydrostatic pressure on individual grains. The aim of the following section is to discuss the influence of nonhydrostatic pressure conditions at the grain level for a polycrystal and the resulting consequences on the constitutive law at the mesoscopic scale.

As a first step, we will focus on the onset of the buckling instability, and thus present in the following a modeling of TATB polycrystal by non linear elasticity in a Lagrangian formulation.

(V /V 0 ) 1/3 I) and the volume variation reads detF = (V /V 0 ). The elastic energy is defined as follows:

with E the Green-Lagrange strain tensor defined by E = 1 2 (F T • F -I) and D(detF ) the volume-variation dependent stiffness tensor. The first Piola-Kirchhoff stress tensor P is obtained through the derivation of E with respect to F :

which after appropriate tensorial derivation leads to:

D(V /V 0 ) can be interpreted as an "isochoric" stiffness tensor, i.e., the projection of the modulus of a certain state of volume/pressure in the reference configuration at constant volume. In order to obtain the expression of D(V /V 0 ), we assume the two following properties:

which implies that D keeps its anisotropy at finite pressure. In order to obtain the expression of D we introduce the pressure-dependent modulus B and its volume-derivative B ′ :

Chapter Summary

• A pseudo-phase transition is observed during MD simulations of TATB under spherical deformation. It cannot be triggered during hydrostatic compression.

• Under spherical deformation, the resulting microstructure is a mix of twinned and pure buckled grains.

• Elastic-Buckling simulations using a element free Galerkin method with non linear elasticity show results in qualitative agreement with MD simulations for the spherical deformation of a single crystal.

• A polycrystalline TATB sample has been simulated through the EB simulation under spherical deformation. The buckling instability is generated at the same compression threshold as MD.

• The average EOS obtained with the polycrystal does not follow any of the EOS computed for the single crystal, indicating the need of a more elaborated modeling of polycrystals with very anisotropic elasticity.

l'échelle du grain d'un polycristal.

Modélisation mésoscopique en élasticité non-linéaire du TATB polycristallin

Une modélisation en élasticité non linéaire en grandes transformations d'un petit polycristal de TATB est proposée. Le but de ce chapitre étant de reproduire le processus de nucléation de l'instabilité de flambage à l'échelle du monocristal, il a été choisi de ne pas prendre en compte le paysage énergétique lié au maclage dans la suite.

Tenseur des modules élastiques en fonction d'une équation d'état

Afin de tenir compte de l'anisotropie élastique du matériau en pression, le tenseur d'élasticité est supposé dépendant de la variation de volume déterminée par une équation d'état. Le tenseur gradient de la transformation est donné par F = (V /V 0 ) 1/3 I et on a la variation de volume detF = V /V 0 . En dérivant le terme d'énergie élastique par rapport au tenseur gradient de la transformation (les détails du calcul sont présentés en section 6.5.1) :

o ù E est l'énergie élastique, E le tenseur de Green-Lagrange et D le tenseur des modules élastiques dépendant de la variation de volume, on obtient une équation différentielle du premier ordre : with angle φ, followed by a rotation θ around the new x-axis, while the third rotation of angle ψ is performed around the new z-axis. However, this formulation leads to a singularity problem since the angular velocities require the differentiation of Euler angles (where ψ and φ become degenerated when θ = 0). In order to overcome this problem, Evans [START_REF] Evans | On the representatation of orientation space[END_REF][START_REF] Evans | Singularity free algorithm for molecular dynamics simulation of rigid polyatomics[END_REF] introduced the use of four parameters as generalized coordinates, through the quaternions. The components of a quaternion in terms of Euler angles are defined as:

with the normalization i q 2 i = 1. The Euler angles based definition of quaternions is very useful to convert to and from the quaternion formal-

with σ PR = Ωh -T a matrix containing informations about size and orientation of the MD simulation cell. The usual set of 3N dynamical variables is then augmented by the nine components of h. The time evolution of these variables is then obtained through the following Lagrangian:

with p the imposed hydrostatic pressure and W a pressure coupling parameter. Minimizing this Lagrangian therefore optimizes both the lattice (through h) and particles locations' to reach the target pressure p.

In this study, a Nosé-Hoover thermostat [START_REF] Nosé | A molecular dynamics method for simulations in the canonical ensemble[END_REF][START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF] has been added in order to control the temperature along with pressure. The idea is to consider a heat bath coupled with the system with an artificial variable ũ, associated with a mass Q > 0 as well as a velocity u. This variable ũ is similar to a time-scaling parameter such that time scale and atomic coordinates are defined as follows:

The Lagrangian for the extended system reads:

where L std is the standard Lagrangian from Equation A.

A.3 Least Square Method for Deformation Analysis

Consider a pair of particles i and j initially separated by a vector ∆X ij from each other, and currently distant by ∆x ij with ∆x ij ≈ F i • ∆X ij . The corresponding error function noted χ 2 ij reads:

Thus, for a number N of neighbors in a sphere of radius r cut around the particle i, the total weighted least squares error is:

where w j is a weight coefficient associated to the j th neighbor varying from 1 to 0 for ∆x ij = 0 to r cut using a cubic spline function. This function is represented in Figure A.1 with abscissa r = ∆x ij /r cut . Here, r cut is the cutoff radius used for the local measure of the deformation and equals two times the cutoff radius of the potential so that a sufficient number of neighbors is taken into accout to perform the measure.

The minimization of χ i with respect to the F components leads to the following relation: