Final thanks go to those people who may not have understood my research but stuck with me and gave me alternative interests outside the lab and office. I will cherish the memorable times spent in their company. Warmest thanks go to my "El-

Résumé

"Etre libre, ce n'est pas seulement se débarrasser de ses chaînes ; c'est vivre d'une façon qui respecte et renforce la liberté des autres."

Nelson Mandela

On peut définir un laser comme étant un appareil qui produit des lumières particulièrement caractérisées par leur cohérence spatiale et temporelle. La cohérence spatiale permet d'obtenir des lumières très étroites et qui peuvent rester étroites sur des longues distances. La cohérence temporelle permet d'obtenir des lumières quasimonochromatiques. Grâce à ces deux propriétés, les lasers sont devenus très utilisés dans des diverses applications couvrant pratiquement tous les domaines. Parmi ces applications, on peut citer le projet laser Mégajoule qui vise à simuler la fusion nucléaire par la technique du confinement inertiel d'une cible de deutérium-tritium (DT) à l'aide des rayonnements lasers très puissants. Dans une telle application, des matériaux peuvent être exposés à des irradiations laser très intenses, ce qui peut causer leur endommagement. Le développement et l'utilisation pratique des lasers de puissance nécessitent des études très approfondies pour bien comprendre les sources et les origines des endommagements causés et éviter les accidents. Ceci a donné naissance à une multitude de travaux expérimentaux qui sont en place depuis les années 60. Malgré les nombreuses avancées faites dans cette direction, la voie expérimentale souffre de plusieurs limitations liées principalement au coût élevé des expériences lasers. Des plus, les difficultés techniques souvent rencontrées empêchent l'obtention d'une information complète de ce qui se passe à l'intérieur du matériau étudié. Ceci a créé le besoin d'un nouvel investissement par la voie numérique pour surmonter les difficultés expérimentales, réduire le coût des études et mieux caractériser l'endommagement laser dans les matériaux. c'est dans ce cadre que s'inscrit ce projet de thèse intitulé "Discrete-continuum coupling method for simulation of laser-induced damage in silica glass".

Le premier objectif de cette thèse est de développer un outil numérique robuste et capable de prédire correctement le comportement mécanique des matériaux soumis à des sollicitations rapides et particulièrement à des sollicitions générées par un rayonnement laser. La plupart des problèmes rencontrés lors des applications lasers sont des problèmes multi-échelles. En effet, une multitude de phénomènes très complexes peuvent interférer dans les zones affectées par le rayonnement laser, où une analyse à très basse échelle est nécessaire pour capter ces phénomènes. Le ix x CONTENTS reste du domaine ne subit que des effets structuraux et une analyse à échelle plus grande suffirait dans cette partie. Une étude multi-échelle pour ce type de problèmes semble donc avantageuse vis à vis du coût et de la précision des résultats. L'analyse des méthodes numériques existantes selon leurs échelles d'application montre qu'elles peuvent être divisées en deux classes : méthodes discrètes basées sur la mécaniques newtonienne (PFD) et méthodes continues basées sur la mécanique des milieux continus (MMC). Les méthodes discrètes semblent bien adaptées pour étudier des problèmes complexes nécessitants une très basse échelle d'analyse tel que le problème de multi-fissuration. Par contre, elles sont très coûteuses en terme CPU. Les méthodes continues sont, en général, moins coûteuses que les méthodes discrètes. Néanmoins, elles sont peu adaptées pour simuler des problèmes complexes impliquant des discontinuités tels que les problèmes rencontrés dans la zone du choc dans les applications lasers. Les avantages de ces deux classes sont largement complémentaires. Alors, dans le but de développer un outil robuste vis à vis des critères de précision et coût de calcul, une méthode multi-échelle de couplage entre une méthode discrète et une méthode continue a été proposée dans cette thèse. Ces deux méthodes (discrete et continue) doivent être rigoureusement choisies pour mieux répondre aux attentes de ce projet "étude d'endommagement des matériaux soumis à des irradiations laser". En effet, un très grand nombre des méthodes discrètes et continues peuvent être trouvées dans la littérature. Chaque méthode a ses avantages et inconvénients. Pour faciliter le choix des approches à coupler, les méthodes numérique ont été classées en différentes catégories. Ceci nous a permis de retenir les deux catégories (discrète et continue) les mieux adaptés à ce projet. Enfin, les méthodes de chaque catégorie ont été comparées pour sélectionnés celles les plus appropriées. Au final, la méthode des éléments discrets (DEM) a été sélectionnée comme méthode discrète et la méthode des éléments naturels contraints (CNEM) à été retenue comme méthode continue. En ce qui concerne les codes de calculs associés, la plate-forme GranOO 1 a été choisie pour effectuer le calcul DEM et la plate-forme NESSY 2 a été choisi pour effectuer le calcul CNEM. Ensuite, un couplage adapté à la dynamique rapide a été développé entre ces deux méthodes (DEM et CNEM). Ce couplage est basé sur la technique d'Arlequin qui procure un vrai contexte physique de couplage de méthodes numériques différentes. Plusieurs paramètres de couplage sont mis en jeu dans cette approche. Dans le but de faciliter le paramétrage de la méthode du couplage DEM-CNEM, et donc son application pratique, une étude paramétrique sur ces différents paramètres a été faite. Suite à cette étude, plusieurs recommandations aidant aux choix de ces paramètres ont été fournies. Enfin pour valider cette méthode de couplage, elle a été appliquée pour simuler quelques tests dynamique de référence. Les résultats obtenus sont très cohérents avec les résultats analytiques et ceux obtenus numériquement en utilisant les deux méthodes, DEM et CNEM, séparément. Sur le plan informatique, cette méthode a été implémentée en trois parties : une partie écrite en C++ pour adapter le code GranOO au couplage, une partie écrite en Python pour adapter le code NESSY au couplage et une partie écrite en Python pour assurer la communication des deux codes GranOO et NESSY. Elle peut CONTENTS xi être utilisé sur une large gamme des problèmes mécaniques muti-échelles. Ceci a fait l'objet d'un premier article qui a été publié dans "Computer Methods in Applied Mechanics and Engineering" [START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF].

Après le développement de l'outil numérique, il est nécessaire de choisir le matériau dont l'endommagement laser sera numériquement investigué. Plusieurs matériaux sont couramment utilisés dans les applications lasers et nécessitent des études supplémentaires de leurs réponses mécaniques aux irradiations lasers. Parmi lesquels, on peut citer le verre de la silice qui, bien qu'il soit le constituant dominant des équipements lasers, est peu étudié numériquement. Ceci est dû à la grande complexité de son comportement mécanique au regard des matériaux classiques tels que les métaux. En effet, il peut subir une déformation permanente même sous pression hydrostatique. En plus, ses propriétés mécaniques dépendent anormalement de la pression et de la température fictive. Dans le but de corroborer les rares études numériques de ce matériau, il a été choisi comme matériau d'étude dans ce travail. Le second objectif de cette thèse est donc de développer un modèle numérique, adapté à la méthode du couplage DEM-CNEM, traduisant le comportement mécanique de ce matériau. La loi de comportement précédemment implémenté dans les deux codes GranOO et NESSY est une loi élastique simple. Cette loi semble suffisante dans les régions CNEM parce que cette méthode est appliquée loin du chargement. Cependant, elle très pauvre pour pouvoir modéliser correctement le comportement mécanique complexe de la silice dans les régions fortement sollicitées qui sont les régions DEM. Pour surmonter cette limitation, deux modèles de comportement adaptés à la DEM ont été proposé dans ce travail. Le premier modèle est basé sur le calcul des contraintes normales dans les joints cohésive reliant les éléments discrets. Ce modèle implique plusieurs paramètres microscopiques à calibrer. Pour faciliter leur calibration, une étude paramétrique a été faite. La validation de ce modèle a été faite, dans un premier temps, en statique pour mieux analyser les éventuels problème en l'absence d'effets dynamiques qui peuvent être complexes. Il a été validé à l'échelle macroscopique via la simulation d'un essai de compression hydrostatique d'une sphère. Les résultats obtenus se comparent favorablement avec des résultats expérimentaux issus de la littérature. À l'échelle microscopique, ce modèle a été validé en reproduisant les essais de micro-indentation de Vickers à faible force d'indentation. Les résultats obtenus sont proches des résultats expérimentaux de la littérature. Ceci a fait l'objet d'un second article qui a été publié dans "Journal of Non-Crystalline Solids" [START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF]. En dynamique, ce modèle a été validé en reproduisant les essais d'impact des plaques en silice. Pour des faible vitesses d'impact, les résultats obtenus se comparent favorablement avec des résultats expérimentaux. Cependant, pour des grandes vitesses d'impact, les résultats numériques, et en particulier l'estimation de la vitesse de l'onde de densification, ne sont pas très satisfaisants. Ceci est dû à l'écart entre la courbe d'Hugoniot numérique et celle obtenue expérimentalement lors des essais d'impact (dont les résultats ont été comparés aux résultats numériques). Cet écart ne peut pas être corrigé avec ce modèle à cause du grand nombre des paramètres à calibrer ce qui rend les simulations plus aléatoires. Pour surmonter ce problème, un autre modèle basé sur la notion du tenseur viriel, équivalent du tenseur des contraintes de Cauchy en MMC, a été proposé. Ce modèle est très flexible et permet de reproduire des courbes d'Hugoniot complexes en utilisant une technique xii CONTENTS de "prédiction-correction". Ces deux modèles ont été implémentés en C++ dans le code GranOO en utilisant l'architecture "Orienté Objet (OO)". En ce qui concerne la modélisation de la rupture, les modèles les plus utilisés dans les méthodes discrètes sont basés sur le calcul d'une contrainte (ou déformation) équivalente dans les joints cohésifs. Si cette quantité dépasse une certaine limite, déterminée par calibration, dans un joint, ce dernier se disparaît. Ces modèles donnent, généralement, des résultats acceptables à l'échelle macroscopique. Cependant, ils n'arrivent pas à reproduire correctement les mécanismes d'endommagement (chemins de fissures) à l'échelle microscopique. Pour corriger ce problème, un autre modèle de rupture basé sur le calcul des contraintes virielles aux éléments discrets a été proposé dans cette thèse. L'application de ce modèle sur quelques problèmes à rupture montre qu'il donne une meilleure représentation des mécanismes d'endommagement à l'échelle microscopique. Ce modèle a été aussi implémenté en C++ dans le code GranOO. Ceci a fait l'objet d'un troisième article qui a été publié dans "Computer Methods in Applied Mechanics and Engineering" [5].

Finalement, en utilisant la méthode de couplage DEM-CNEM et les différents modèles de comportement et de rupture, l'endommagement induit par un rayonnement laser de la silice a été étudié via la simulation de l'essai choc laser sur la silice. Par manque des résultats expérimentaux, seuls des résultats qualitatifs ont été présentés dans cette étude. Qualitativement, les résultats obtenus sont très cohérents avec les observations expérimentales et la théorie du choc.

Acronyms

Motivation and objective research

The term "laser", acronym for Light Amplification by Stimulated Emission of Radiation, is a device that emits light through a process of optical amplification based on the stimulated emission of electromagnetic radiation. Lasers differ from other sources of light, because they emit spatial and temporal coherent lights. Their spatial coherence allows them to be focused to a tight spot and keep a laser beam collimated over long distances [START_REF] Haken | Laser theory[END_REF]. Their temporal coherence allows them to have a very narrow spectrum (i.e. virtually a single light color is emitted) and to emit light pulses that only last a few femtoseconds [START_REF] Haken | Laser theory[END_REF]. Because of these excellent coherence properties, lasers have become very used in several low-power laser applications such as the low-level laser therapy (LLLT) in medicine. Great success and increased variety of laser devices have enabled an extension for the use of high-power laser radiations in scientific, industrial, and military applications. Of particular applications are the Laser Mégajoule (LMJ) at the Atomic Energy and Alternative Energies Commission (CEA) in France 1 and the National Ignition Facility (NIF) [START_REF] Brereton | Overview of the National Ignition Facility[END_REF][START_REF] Lindl | The physics basis for ignition using indirect-drive targets on the National Ignition Facility[END_REF] at the Lawrence Livermore National Laboratory (LLNL) in the United States of America, which are designed to achieve the inertial confinement fusion (ICF) of a deuterium-tritium (DT) target. In these applications, 192 laser beams will produce 1.8 MJ of energy aimed at a small target containing a spherical DT pellet (Fig. 1.1). This process produces a large amount of heat and generates local temperatures comparable to the sun and other stars. Such temperatures enable the fusion of the deuterium and tritium nuclei and releasing more energy than was put in by laser. Figure 1.2 summarizes the different steps of the ICF. In these applications, materials would be subjected to laser radiations, which can lead to serious damages. Indeed, when a material is irradiated by a laser beam, compressive stresses are generated by means of the recoil pressure due to the rapid expansion of the plasma plume resulting from the interaction between the first atomic layers and the laser beam (Fig. 1.3). The propagation of these compressive stresses within the material can damage it. Development and practical use of high-energy lasers therefore require in-depth studies of the origins and sources of the laser-induced damage in materials. Such research efforts have been in place since the sixties and have resulted in development and commercial manufacturing of variety of specialized materials that can be used safely in high-power laser systems. Numerous experimental techniques and processes have been proposed to study the laser-induced damage in materials. Among them, one can cite the Laser Shock Processing (LSP) which is based on the focusing of a laser beam on a specimen and the detection of some irreversible damage in it. Despite the current state of the experimental studies, the laser-induced damage in materials is always one of the major concerns encountered in high-power laser applications. In effect, the experimental works are very costly. Besides, technical difficulties presently prevent the effectiveness of experimental works in terms of providing insightful and complete Figure 1.3: Laser-matter interaction information that cannot be directly measured, observed, or is difficult to acquire via experimental means. Thus, a next generation of research efforts is required and the numerical support is necessary to overcome the experimental difficulties, to reduce the cost of studies and to better characterize the laser-induced damage in materials.

The first objective of this dissertation is to develop a robust numerical tool able to predict correctly the mechanical response as well as the damage mechanisms of materials subjected to dynamic stresses and particularly to stresses generated by laser radiations.

After developing the numerical tool, it is necessary to choose the material in which the laser-induced damage will be numerically studied in this work. Several interesting materials are routinely used in different laser applications and require additional investigations in their responses to laser radiations. Of particular interest is the silica glass which is the dominant constitutional material of the optical equipments in laser devices. This material is known to exhibit anomalous behavior in its thermal and mechanical properties [START_REF] Brückner | Properties and structure of vitreous silica. I[END_REF][START_REF] Brückner | Properties and structure of vitreous silica[END_REF]. Furthermore, certain properties of this glass such as Young's modulus, shear modulus and density show anomalous dependence on the fictive temperature. Because of its complex mechanical behavior, numerical study of this material remains a central issue for several researchers. These reasons have made silica glass an attractive material to be numerically studied. Therefore, this material is chosen in this work to numerically investigate its mechanical response to high-power laser radiations. The second objective of this dissertation is to model the silica glass mechanical behavior using the numerical tool which will be developed in this work.

Finally, the mechanical part of the Laser Shock Processing (LSP) on silica glass using high-power laser will be reproduced numerically to study the laser-induced damage in this material.

Methodology of the work

As a first part of this work, a new numerical tool is developed, providing an alternative way of scientific investigation, instead of carrying out expensive, timeconsuming or even dangerous experiments in laboratories.

Because of the complexity of the spawned phenomena in laser processes [START_REF] Kruusing | Underwater and water-assisted laser processing: Part 1-general features, steam cleaning and shock processing[END_REF][START_REF] Kruusing | Underwater and water-assisted laser processing: Part 2-Etching, cutting and rarely used methods[END_REF][START_REF] Berthe | State-of-the-art laser adhesion test (LASAT)[END_REF][START_REF] Fabbro | Physics and applications of laser-shock processing[END_REF][START_REF] Cottet | Spallation studies in aluminum targets using shock waves induced by laser irradiation at various pulse durations[END_REF][START_REF] Boustie | Experimental and numerical study of laser induced spallation into aluminum and copper targets[END_REF][START_REF] Tollier | Study of the laser-driven spallation process by the velocity interferometer system for any reflector interferometry technique. I. Laser-shock characterization[END_REF][START_REF] Tollier | Study of the laser-driven spallation process by the VISAR interferometry technique. II. Experiment and simulation of the spallation process[END_REF], the continuum methods, such as the Finite Element Method (FEM) [START_REF] Zienkiewicz | Finite Element Method: Its Basis & Fundamentals[END_REF][START_REF] Zienkiewicz | The Finite Element Method for Solid and Structural Mechanics[END_REF][START_REF] Zienkiewicz | The Finite Element Method for Fluid Dynamics[END_REF] and the Constrained Natural Element Method (CNEM) [START_REF] Chinesta | Natural element method for the simulation of structures and processes[END_REF],

are not well suited to simulate these processes. In effect, this class of methods is based on the continuum mechanics, and additional handling must be performed to take into account the prospective micro-discontinuities generated by the cracking mechanisms. The discrete methods, such as Molecular Dynamics (MD) [START_REF] Muralidharan | Molecular dynamics studies of brittle fracture in vitreous silica: Review and recent pregress[END_REF] and Discrete Element Method (DEM) [4], seem more adapted to treat these problems. Because, they naturally take into account the discontinuities. The major drawback of these methods is that they are very time-consuming (CPU-wise) and the computation time can become crippling, especially in the case of large studied domains. However, in most situations and particularly in the case of LSP processing, the regions requiring fine scale analyses by discrete methods are small with respect to the full studied domain. Thus, the use of a specific multi-scale method to treat phenomena at each scale seems to be advantageous. This is why a discrete-continuum coupling method, adapted for dynamic simulations, is proposed in this dissertation. First, a bibliographic review of the most used numerical methods is performed to choose the appropriate continuum and discrete methods to be coupled. Following this review, the Constrained Natural Element Method (CNEM) [START_REF] Chinesta | Natural element method for the simulation of structures and processes[END_REF] is chosen as the continuum method, whereas, the Discrete Element Method (DEM) [4] is chosen as the discrete method. Second, a discrete-continuum coupling approach, adapted for dynamic simulations, between these methods (DEM and CNEM) is developed. This coupling method is based on the Arlequin approach [START_REF] Dhia | Analyse mathématique de la méthode Arlequin mixte[END_REF][START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF] and involves several coupling parameters. To simplify the application of the new discrete-continuum coupling method on complex applications, a parametric study of the different coupling parameters is performed. At the end of this study, several practical recommendations that simplify the coupling parameters setting in practice are given. The validation of this approach is performed using several dynamic reference tests. The developed coupling method could be applied to a wide range of dynamic problems in which the fine scale effects are localized in a small portion compared to the whole studied domain. This part has been the subject of a first paper [START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF] published in Computer Methods in Applied Mechanics and Engineering.

The DEM-CNEM coupling method will be used to simulate the LSP processing on silica glass. The shocked zone, where complex phenomena can take place, will be modeled by DEM, whereas the remaining estate will be modeled by CNEM. A question that arises here is how to model the complex dynamic behavior of silica glass using this approach. Because the CNEM is applied far from the high loaded regions, an elastic constitutive law is used in the associated regions. However, a complex constitutive law taking into account the different domains of the silica glass mechanical behavior (elasticity, densification and saturation) must be used in the DEM regions. As a second part of this work, a new beam-based complex model of silica glass mechanical behavior is developed for the DEM. This model is first validated in statics to better analyze the potential problems in the absence of complex dynamic effects. Spherical samples subjected to hydrostatic pressures are used to validate this model at macroscopic scale. Then, simulation of Vickers micro-indentation is performed for the microscopic validation. This work has been the subject of a second paper [START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF] published in Journal of Non-Crystalline Solids. Once statically validated, this model is validated in dynamics by simulation of high velocity impacts of plates. While this model provides relatively good results, its application remains difficult. Because it involves several microscopic parameters which must be determined by calibration. To overcome this limitation, another model based on the computation of virial stress [START_REF] Mclellan | Virial theorem generalized[END_REF][START_REF] Zhou | A new look at the atomic level virial stress: on continuum-molecular system equivalence[END_REF], equivalent to Cauchy stress in continuum mechanics, at the discrete elements is developed. Reproducing the high velocity impact simulations using this new model, good results are obtained, compared to those obtained using the first model.

To model fracture in discrete methods, the approaches most commonly used are based on the computation of bond's strain [START_REF] Carmona | Fragmentation processes in impact of spheres[END_REF] or stress [START_REF] Potyondy | Abonded-particle model for rock[END_REF]. The bond breaks if the associated strain (or stress) exceeds a given value determined by calibration tests. At the macroscopic scale, the fracture results of these criteria are moderately satisfactory [3]. However, at microscopic scales, these criteria cannot reproduce correctly the cracking mechanisms [3]. In this work, a new fracture model adapted for the Discrete Element Method is developed. This model is based on the computation of the virial stress [START_REF] Mclellan | Virial theorem generalized[END_REF][START_REF] Zhou | A new look at the atomic level virial stress: on continuum-molecular system equivalence[END_REF] at each discrete element. As will be seen, this model reproduces well the cracking mechanisms at both macroscopic and microscopic scales, and corrects the problems encountered using the existent models. This part has been performed in collaboration with D. André and has been the subject of a third paper [5] published in Computer Methods in Applied Mechanics and Engineering.

Finally, all these developments are used to simulate the Laser Shock Processing on silica glass. Due to lack of enough experimental results, only qualitative results of this simulation are given. Qualitatively, numerical results compare well to experimental observations reported in the literature.

Figure 1.4 summarizes the objectives of this dissertation. As seen in this figure, the main request of this dissertation, Discrete-continuum coupling method for simulation of laser-induced damage in silica glass, is divided into four steps as follows:

• development of a robust numerical tool adapted for dynamic simulations,

• development of a constitutive model taking into account the different domains of the silica glass mechanical behavior,

• development of a fracture model able to reproduce correctly the microscopic cracking mechanisms,

• simulation of Laser Shock Processing on silica glass to study the mechanical damage induced by such a loading.

Dissertation structure

Following this introduction, this dissertation will be laid as follows:

• Chapter two will introduce some background knowledge of the different disciplines that interact to accomplish the objectives of this dissertation: first, the laser-matter interaction theory will be briefly recalled; then, the silica glass mechanical behavior will be detailed; finally, the main specificities of the numerical simulation and the different classes of the numerical methods will be reviewed.

• Chapter three will review the principle and the main specificities of the most commonly used numerical methods to simulate the mechanical behavior of materials. The focus of this chapter is to choose the appropriate discrete and continuum methods to be coupled in this work.

• Chapter four will detail how the retained discrete and continuum methods are coupled. Then, a parametric study of the different coupling parameters will be made to draw recommendations simplifying the choice of these parameters in practice and so the application of the coupling method on complex problems. Finally this method will be validated using several dynamic reference tests: dynamic tension, bending, torsion and impact tests on 3D beams.

• Chapter five will discuss how the mechanical behavior of silica glass is modeled. Two models are proposed in this chapter. The first is based on the normal stress in the cohesive beam bonds. This model is first validated in statics to simplify the analyses of the potential difficulties, and then in dynamics by simulation of high velocity impacts of silica glass plates. The second model is based on the "prediction-correction" technique, using the virial stress at the discrete elements. As will be seen later, this model is easier to be applied and gives good results compared with the first model. Subsequently, it will review how the brittle fracture of silica glass is modeled.

• Chapter six will study mechanically the laser-induced damage in silica glass.

To do this, the Laser Shock Processing on silica glass is reproduced numerically using the new discrete-continuum coupling method as well as the models of silica glass behavior and brittle fracture developed in this work. Due to the lack of sufficient experimental results, only a qualitative study is given in this chapter.

• Finally, this dissertation will be closed by several conclusions and outlooks.

Laser-matter interaction theory

The interaction mechanisms between matter and electromagnetic waves depend on the waves frequency range. The electromagnetic spectrum (Fig. The interaction between a laser beam and matter is based on the fast photon energy change into thermal and/or kinetic energy in the first atomic layers of matter. When a material is irradiated with lasers, the laser energy will first be transformed into electronic excitation energy and then transferred to lattices of the material through collisions between the electrons and lattices. The deposition of laser energy will produce a series of effects, such as temperature rise, gasification and ionization. The physical processes of laser-matter interaction (Fig. 2.2) depend mainly on the laser intensity

I L = E L τ L S L
, where E L is the pulse energy, τ L is the Full Duration at Half Maximum (FDHM) and S L is the energy deposition surface. The present work will focus on laser radiations with laser density above 10 12 W cm -2 , leading to plasma generation and ablation of matter.

Generation of shock wave by laser ablation

Materials subjected to laser irradiation (Fig. 2.3) will absorb the incident laser energy [START_REF] Han | Interaction Between Pulsed Laser and Materials[END_REF]. In the case of normal incidence, the coefficient of energy absorption A is given by Equation (2.1):

A = 4 n L (n L + 1) 2 + k 2 L (2.1)
where n L and k L are respectively the real and imaginary parts of the refraction complex index ( ñL = n L + ik L ). Otherwise, the influence of the incidence angle and the polarization of the laser radiation must be taken into account in the computation of A. The laser energy is firstly absorbed by the free electrons in a small depth l skin determined by Equation (2.2):

l skin = λ 2 π k L (2.2)
with λ is the wavelength of the incident laser wave. Then, the absorbed energy is The energy absorption leads, inter alia, to material heating. In the case of normal incidence, the temperature raise can be expressed in terms of depth x and laser irradiation time t as follows (2.3):

∆T(x, t) = 2 (1 -R) α I L t π k ρ C ierfc   x 2 k t ρ C   (2.3)
where R is the reflectivity, α is the absorptivity, I L is the laser intensity, k is the thermal conductivity, C is the specific heat, ρ is the density of the irradiated material and "ierfc" is the inverse complementary error function. According to τ L , two types of laser pulses can be distinguished: short pulses if τ L < τ ei and long pulses if τ L >= τ ei . In the case of a short pulse, the energy absorption and the material heating are decoupled (ions are not heated during the laser irradiation). The absorbed energy is transfered to ions only from the free electrons by collision. If the ions energy exceeds the ionization threshold, the heated zone is transformed to a plasma. Whereas, in the case of a long pulse, these phenomena are coupled and the laser radiations interact directly with the plasma. This interaction consists of :

• the temperature of the irradiated material will rise and the deposited energy diffuses inside the material to a certain depth (Fig.

2.4a);

• the temperature of melted material sharply increases to over the boiling point due to the heavy deposition of laser energy; nevertheless, the boiling does not start and the liquid is super-heated because of the absence of nucleation (Fig.

2.4b);

• the disturbance will bring about nucleation and the super-heated liquid thickens as the size and the number of bubbles grow (Fig.

2.4c);

• the startling boiling will arise once the size of bubbles is sufficiently large and afterwards the super-heated liquid and particles will be ejected (Fig.

2.4d).

This process is known as material laser ablation.

As a consequence of the material laser ablation, high pressures will raise and propagate in the material. The maximum value of the laser-induced pressure P max L can be determined using experiment-based abaci [START_REF] Lescoute | Étude de la fragmentation dynamique de métaux sous choc laser[END_REF]. In nanosecond regime, P max L can be approximated using the empirical formula of Phipps [START_REF] Phipps | Laser ablation powered mini-thruster[END_REF] as follows:

P max L = 0.622 10 11 A 7/16 Z 9/16 (α I L ) 3/4 (λ √ τ L ) 1/4 (2.4)
where A and Z are respectively the mass and atomic number of the irradiated material. However, it is very difficult to determine experimentally the temporal profile of the laser-induced pressure in materials which is necessary to model any laser application. This one can be estimated numerically by solving the equations of lasermatter interaction. Several software specialized in such a computation can be found in the literature, one can cite for example ESTHER [START_REF] Colombier | Optimized energy coupling at ultrafast laser-irradiated metal surfaces by tailoring intensity envelopes: Consequences for material removal from Al samples[END_REF] and CHIC [START_REF] Maire | A cell-centred arbitrary Lagrangian-Eulerian (ALE) method[END_REF]. The former is developed in CEA/DIF Île-de-France and allows 1D simulations of interaction between light and matter, based on the Helmhotz equations [START_REF] Khan | Wavelength-Dependent Modifications in Helmholtz Optics[END_REF]. Therefore, it can be used to perform computation in the femtosecond regime. The latter is a bidirectional (2D) code developed in the CELIA laboratory. The equations of light-matter interaction implemented in this code are based on the theory of collision-induced absorption [START_REF] Frommhold | Collision-induced Absorption in Gases[END_REF] which limits its application to the nanosecond regime. After determining the temporal profile of the laser-induced pressure, it will be possible to study the mechanical effects of this loading which, under certain conditions, can lead to a shock wave. This part of laser-matter interaction is not within the scope of this dissertation which aims to simulate the mechanical behavior of silica glass under Laser Shock Processing (LSP). The temporal pressure profile is assumed to be given and will directly be applied on the shock surface (irradiated surface). Also, it is assumed that the thermal effects are relatively slow compared to the propagation of the laserinduced pressures within the material. Therefore, the thermal aspect of the LSP processing will be neglected in the present work.

Shock wave propagation in materials

The laser-induced pressure will make the irradiated material under compressive stresses which will propagate in all directions [START_REF] Tani | Warm Laser Shock Peening: New developments and process optimization[END_REF][START_REF] Lescoute | Étude de la fragmentation dynamique de métaux sous choc laser[END_REF][START_REF] Fabbro | Physics and applications of laser-shock processing[END_REF]. When the amplitude of the stress waves greatly exceeds the dynamic strength of the material, the shear stresses can be neglected in comparison with the hydrostatic component of the stress tensor. One therefore can consider a high pressure state traveling into the material which can be assumed, in a first approximation, to have no shear resistance (i.e. the shear modulus is zero, µ = 0). Under this assumption, its state is completely characterized by three thermodynamic parameters: pressure P, density ρ (or specific volume v = ρ -1 ) and internal energy E (or temperature T). The propagation of a high pressure state into a material can lead to the formation of a shock wave. To simplify the understanding of this concept, an ideal gas will be first considered.

For an ideal gas, the associated equation of state in the case of an isentropic process can be written:

PV γ = const (2.5)
where V is the volume occupied by the gas, γ = C p C v is the adiabatic index, C p and C v are the specific heat at constant pressure and the specific heat at constant volume, respectively. Derivation of Equation (2.5) leads to:

∂P ∂V = -γ P V (2.6) Equation (2.6) shows that the compressibility K = - 1 V
∂V ∂P s of an ideal gas decreases with pressure P, since P V increases with P. This is explained physically by the fact that the external electron shells of adjacent atoms penetrate and repel each other. Therefore, the velocity of a disturbance C in an ideal gas, which is equivalent In other words, the high amplitude isentropic disturbances travel faster than the low amplitude ones. This is the sine qua non requirement for shock wave formation. Therefore, a compression front will steepen up as it travels through the material because the higher amplitude regions of the front travel faster than the lower amplitude regions. This leads to a shock wave which is defined simply as a discontinuity in the thermodynamic parameters (P, ρ (or v) and E (or T)). Conversely, a release front will spread out during its propagation through the material (Fig.

2.5).

For a solid, it is necessary to differentiate between the deviatoric and hydrostatic components of stress. When the former is negligible, it becomes possible to apply the treatment for fluids without remorse.

Hydrodynamic treatment

Shock waves are characterized by a steep front and require a state of uniaxial strain (no considerable lateral flow of material) which allows the buildup of the hydrostatic component of stress to high levels. When the last (hydrostatic component) far exceeds the dynamic flow strength, the material behaves as a fluid. Therefore, the Rankine-Hugoniot conservation equations for fluids [START_REF] Rankine | On the Thermodynamic Theory of Waves of Finite Longitudinal Disturbance[END_REF][START_REF] Rankine | Sur la propagation du mouvement dans les corps et spécialement dans les gaz parfaits, 1e Partie[END_REF][START_REF] Rankine | Sur la propagation du mouvement dans les corps et spécialement dans les gaz parfaits, 2e Partie[END_REF] can be applied to calculate the shock wave parameters. This is valid taking into account several assumptions:

• the shock is a discontinuous surface and has no apparent thickness;

• the material behaves as a fluid (µ = 0); the theory is, therefore, restricted to high pressures;

• body forces (such as gravitational) and heat conduction at the shock front are negligible;

• there is no elastoplastic behavior;

• material does not undergo phase transformation.

Taking into account these assumptions, the shock wave equations can easily be obtained by considering a small region Ω s with cross section A immediately ahead of and behind the shock front (Fig.

2.6).

Ahead of the front (initial state), the pressure is P 0 , the density is ρ 0 and the energy is E 0 ; behind it, they are P, ρ and E, respectively. The velocity of the particles Figure 2.6: Illustration of a shock front ahead of and behind the shock front, which is moving at a velocity of V s (shock velocity), are (V p ) P 0 = V 0 and V p , respectively. The apparent velocity of the shock front is (V s -V 0 ), because it is moving into a region of particle velocity V p = V 0 . At the same time, the material leaving the shock front is moving at a velocity V s -V p . Using this, the equations of mass, momentum and energy conservation can be derived:

Mass conservation:

This principle states that the mass entering Ω s equals the mass leaving it, hence:

A ρ 0 (V s -V 0 ) dt = A ρ (V s -V p ) dt (2.7)
which can be written simply as:

ρ 0 (V s -V 0 ) = ρ (V s -V p ) (2.8)

Momentum conservation:

The momentum change in Ω s must be equal to the impulse given to it:

ρA(V s -V p )dt mass V p Particle velocity -ρ 0 A(V s -V 0 )dt mass V 0 Particle velocity = (P -P 0 ) A dt Impulse (2.9)
Taking into account (2.8), Equation (2.9) can be written:

ρ 0 (V s -V 0 )(V p -V 0 ) = P -P 0 (2.10)

Energy conservation:

The conservation of energy postulates that the work of the external forces (pressure in this case) on Ω s (∆W) equals the difference in the total energy (kinetic and internal) in Ω s (∆E tot ):

∆W = ∆E tot (2.11) with: ∆W = PA Force V p dt Distance -P 0 A Force V 0 dt Distance (2.12)
and

∆E tot = 1 2 ρ A (V s -V p ) dt V 2 p + E A ρ (V s -V p ) dt - 1 2 [ρ 0 A (V s -V 0 ) dt] V 2 0 -E 0 A ρ 0 (V s -V 0 )dt (2.13)
Taking into account (2.8), Equation (2.14) can be expressed:

P V p -V 0 = ρ 0 (V s -V 0 ) E -E 0 + 1 2 V 2 p -V 2 0 (2.14)
By summation of (2.10) and (2.14) and introducing (2.8), one can obtain the classic form of the energy conservation:

E -E 0 = 1 2 (P -P 0 ) (v 0 -v) with v = 1 ρ (2.15)
Equations (2.8), (2.10) and (2.15) are known as the Rankine-Hugoniot equations. These equations involve five parameters which are the pressure P, particle velocity V p , shock velocity V s , density ρ and energy E. Hence, an additional equation is required to determine these parameters as a function of one of them. This fourth equation which can be conveniently expressed as the relationship between the particle and shock velocities is called the equation of state and has to be experimentally determined. An empirical polynomial equation with parameters C 0 , S 1 , S 2 , etc, is proposed in the literature to describe this relationship (2.16).

V s = C 0 + S 1 V p + S 2 V 2 p + ... (2.16)
With C 0 is the hydrodynamic sound velocity in the material at zero pressure, expressed in terms of longitudinal C l and transverse C t sound velocities as

C 0 = C 2 l - 4 3 C 2 t ,
and S 1 and S 2 are empirical parameters. This equation (2.16) is often called the equation of state (EOS) of a material. For most materials in the absence of phase transformation, this equation reduces to a linear relationship (S 2 = 0) as shown in Figure 2.7. Several works giving the values of the EOS parameters for different materials can be found in the literature [START_REF] Meyers | Shock waves and high-strain-rate phenomena in metals : concepts and applications[END_REF][START_REF] Walsh | Equation of State of Metals from Shock Wave Measurements[END_REF]. Applying this equation with the conservation equations, other relationships between these parameters can be determined, e.g. P -V s , P -V p , Pρ and E -V s relationships.

Representative curves for shock waves

Hugoniot curve: The relationship between P and ρ (or v = ρ -1 ) is usually known as Hugoniot equation. This equation is defined as the locus of all shocked states in the (P, ρ) or (P, v) frame and essentially describes the material properties (Fig. 2

.8a).

The straight line in Figure 2.8a relaying the (P 0 , v 0 ) and (P 1 , v 1 ) states is the Rayleigh line and refers to the shock state at P 1 . It is very important to realize that when pressure is increased in a shock front, it does not follow the Pv curve. Rather, it jumps discontinuously from P 0 to P 1 . The slope of this line is proportional to the square of the shock velocity (V s ). (Taken from [START_REF] Meyers | Dynamic behavior of materials[END_REF])

Shock polar:

The shock polar (Fig. 2.8b) is defined as the locus of all shocked states in a material expressed in the (P, V p ) frame. In this plot, the slope of the Rayleigh line is proportional to the shock velocity (V s ) and called the shock impedance: Z = ρ 0 V s . In the region of pressures corresponding to low particle velocities in comparison with the shock velocity (V p ≪ V s ), the shock impedance is reduced to the acoustic impedance Z = ρ 0 C 0 . In this case, the shock polar can be approximated by a straight line.

Permanent deformation and shock wave

The above description of shock waves are given such as the shock front is considered as a hydrostatic pressure P. Thus, in a Cartesian coordinate system, the associated stress tensor can be written:

σ =   P 0 0 0 P 0 0 0 P   (2.17)
This is particularly true in the case of pressures widely exceeding the dynamic material strength. In the region of low shock pressures, the irreversible effects which raise during the shock wave propagation become significant. These effects are illustrated, hereafter, using a simplified elastoplastic model based on the Tresca criterion. When a monodimensional shock wave is traveling into an homogeneous and isotropic solid, a state of uniaxial strain is established. Under these conditions, the (a) Hugoniot: Pressure (P) vs. specific volume (v) (b) Shock polar: Pressure (P) vs. particle velocity (V p ) Figure 2.8: Representative plots for shock waves deviatoric component of the stress tensor is not zero. Assuming that the shock wave is traveling with respect to xaxis, the stain and stress tensors are written, respectively:

ε =   ε x 0 0 0 0 0 0 0 0   (2.18) σ = σ s + σ d =   Kε x 0 0 0 Kε x 0 0 0 Kε x   +   4 3 Gε x 0 0 0 -2 3 Gε x 0 0 0 -2 3 Gε x   (2.19)
where K and G are respectively bulk and shear moduli of the material. Under the dynamic elastic limit, called the Hugoniot Elastic Limit HEL, the solid behaves in a perfectly elastic manner. Beyond this value, it begins to exhibit signs of permanent deformation (plasticity). In terms of material strength and using the Tresca formalism, a material deforms plastically when the shear stress τ max equals the yield shear stress τ 0 :

τ max = 1 2 max i,j σ i -σ j = τ 0 (2.20)
where σ i , i = 1, 2, 3 are the principle stresses. Using (2.20), the relationship between HEL and τ max can be derived: or :

P = 1 3 (σ 1 + σ 2 + σ 3 ) = 1 3 σ x + 2σ y = σ x - 2 3 σ x -σ y
τ max = 1 2 σ x -σ y (2.21)
thus,

σ x = P + 4 3 τ max =⇒ HEL = P HEL + 4 3 τ 0 (2.22)
When the HEL is reached, the Hugoniot curve shows a change in slope (Fig. 2.9). Furthermore, the elastic portion of the wave is separated from the plastic wave (Fig. 2.9). One can distinguish two shock front : elastic shock front traveling at the longitudinal sound velocity C l and plastic shock front traveling at the shock velocity V s . The consequences of this phenomenon will be discussed in the next subsection ( §2.1.2.4). At sufficiently high pressures, the plastic shock front overcomes the elastic one, since V s increases with P. This justifies the application of the shock waves hydrodynamic theory in this region of pressures.

Phenomena spoiling the propagation of shock waves

Hydrodynamic damping: As shown above, for most materials, the concavity of the Hugoniot curve is positive. Therefore, the compression front will steepen up as it travels through the material because the higher amplitude regions of the front travel faster than the lower amplitude regions. In contrary, the release front will spread out (the release rate decreases) as it travels through the material. The head of the release part (at the maximum pressure) is traveling at a velocity (V p + C) P f which is greater than the shock velocity V s , where V p and C are the particle velocity and the longitudinal sound velocity at the pressure P f . As the wave progresses, the release part of the wave overtakes the front. This will reduce the pulse duration to zero. After it is zero, the peak pressure starts to decrease and so the shock velocity (Fig.

2.10)

. This is generally referred to as "hydrodynamic damping".

Elastoplastic damping:

In addition to the hydrodynamic damping, the shock waves in solids can also be attenuated due to the difference between the velocity of the elastic and plastic fronts. As seen previously, beyond the HEL, the shock wave is split into two parts: elastic shock front and plastic shock front traveling respectively at velocities C l and V s . During rarefaction (release), the elastic part of the release wave follows the shock front at velocity C l which is, in the region of low pressures, greater than V s . In this case, it can overtake the shock front and reduce the peak pressure (Fig. 2.10). This mechanism of attenuation is classically known as "elastoplatic damping".

Shock wave transmission and reflexion

During its propagation, the shock wave can migrate between media, e.g. when reaching the free surface (material/air interface) in the case of LSP processing, which can result in changes in thermodynamic and kinetic parameters. This depends mainly on the shock impedance of the considered media (Z = ρ 0 V s ). The shock impedance of materials is all the higher as they have high sonic velocity and density. The best way to study the wave transfer between two media is by means of the impedance matching technique. The continuity at the materials interface dictates that the particle velocity and pressure must be the same in both materials. Considering the case of a shock wave of pressure P 1 propagating between two media A and B having impedances Z A and Z B , two cases can be distinguished:

1. Z A < Z B : At the interface, the incident pressure P 1 will change to ensure the dynamic equilibrium of this region which can graphically be determined using the shock polars of A and B. After plotting these curves, the reflected curve (AR) can be plotted such as it passes through the shock pressure P 1 (Fig.

2.11)

. The last curve (AR) intersects the curve B at P 2 . This is the pressure in medium B. As the shock front reaches the interface, the pressure rises to P 2 and propagates into both A and B. The same reasoning applies to determine the particle velocity in both A and B at the interface V P 2 . The bottom view of Figure 2.11 presents the associated sequences of pressure profiles.

Z A > Z B :

A similar reasoning as in the previous case can be applied with reversing the curves A and B (Fig. 2.12). In this case, P 2 is lower than P 1 (P 2 < P 1 ) which produces a release pulse to be sent through the medium A. It should be noticed that the release portion of the wave is not a shock wave but has a slope that decreases with increasing propagation velocity.

The second case (Z A > Z B ) is frequently encountered in the LSP experiments, which explains the spalling phenomenon often observed at the rear surface of the irradiated Figure 2.11: Transmission of shock wave from material A with low impedance to material B with high impedance material. Indeed, when the laser-induced shock wave reaches this surface (if not coated, material/air interface), a release pulse is produced and sent through the irradiated material. This pulse will travel freely until it encounters the release portion of the primary pulse. At this time (t 4 in Figure 2.12), a tensile pulse will form which will propagate in both directions. If this tensile pulse has sufficient amplitude, a spall can be formed. This phenomenon will be studied in Chapter 6.

Laser-induced damage in materials

When a material is irradiated by a laser beam, there are three main classes of mechanisms which give rise to laser-induced damage [START_REF] Wood | Laser-Induced Damage of Optical Materials[END_REF]:

• Thermal processes: they arise from absorption of the laser energy in the material, and in general they occur for continuous wave operation, long pulse lengths and high-pulse-repetition-frequency pulse trains.

• Dielectric processes: they arise when the electric field density is high enough to strip electrons from the lattice. They take place when the pulse lengths are short enough for avalanche ionization to take place and when the thermal absorption is low enough for the avalanche threshold to be below the thermal threshold.

Figure 2.12: Transmission of shock wave from material A with high impedance to material B with low impedance

• Mechanical processes: by propagation of a non-sustained shock wave into a solid, it is possible to produce strong levels of tensile leading to the damage of the material [6]. The laser driven shock wave propagates into the material and reverberates on the back surface (free surface) into a release wave which interacts with the incoming release wave due to the fast laser unloading. The crossing of these two release waves results in tensile state inside the material. The distance between the tensile zone and the back surface is proportional to the pulse duration. If the amplitude and the duration of the tensile state are sufficient, there can be fracture of the considered material and formation of a spall (fig. 2.13a). The incident wave continues to reflect on the new free surfaces (generated by spalling) and the fracture threshold σ f can be reached another time. Therefore, second spall can be generated (Fig. 2.13b). The spalling process will continue as long as the incident wave is sufficiently intense. This is known as multiple spalling (or multiple scabbing). Each spall is generated from a lower compression state, and thus with a lower initial velocity. The multipe scabbing is clearly evidenced in Figure 2.14 taken from the experimental works of Rinehart [START_REF] Rinehart | Scabbing of Metals under Explosive Attack: Multiple Scabbing[END_REF].

The first two processes are due to laser energy absorption (first step of lasermatter interaction) and generally occur in a small depth from the front face of the irradiated material. These processes are not within the scope of this dissertation which focuses on the laser-induced damage from a mechanical point of view. Hence, this work focuses only on the third process which will be studied numerically in Chapter 6.

Mechanical behavior of silica glass

According to the dependence of various properties such as density, hardness, refractive index, Young's and shear moduli with the pressure and fictive temperature (defined in a simple manner as the temperature at which the glass would be in thermal equilibrium), glasses can be classified under three headings: normal, intermediate and anomalous. Normal glasses, such as soda-lime silicate glass [START_REF] Houérou | Surface damage of soda-lime-silica glasses: indentation scratch behavior[END_REF][START_REF] Rouxel | The brittle to ductile transition in a soda-limesilica glass[END_REF], exhibit properties decreasing with increasing fictive temperature and pressure. Anomalous glasses, such as germania, exhibit the opposite trend in these properties. Interme-

Young's modulus Poisson's ratio Density

E = 72 GPa ν = 0.17 ρ = 2200 Kg m -3
Table 2.1: Elastic properties of silica glass diate between normal and anomalous glasses are expected to have some properties that are independent of these parameters (i.e. pressure and fictive temperature). Silica glass, which is a typical amorphous material, belongs to the third heading: anomalous glasses. This material presents an anomalous behavior in its thermal and mechanical properties [START_REF] Brückner | Properties and structure of vitreous silica. I[END_REF][START_REF] Brückner | Properties and structure of vitreous silica[END_REF] which generally increase with pressure and fictive temperature. This section will detail the main specificities of the silica glass mechanical behavior. To simplify the comprehension of this behavior, the silica glass response under static hydrostatic compression will be first studied. And after, the particularity of this behavior under shock compression will be seen.

Silica glass response under static hydrostatic compression 2.2.1.1 Elasticity

In the region of ordinary pressures (up to approximately 8 GPa), the silica glass behaves in a perfectly elastic manner. As mentioned in several experimental works [31,[START_REF] Bridgman | Rough compressions of 177 substances To 40,000 KG/CM2[END_REF][START_REF] Kondo | Nonlinear pressure dependence of the elastic moduli of fused quartz up to 3 GPa[END_REF], for relatively low pressures (under 3 GPa), this material shows a nonlinear elastic behavior which can be explained by the pressure dependence of the bulk and shear moduli, noted respectively K and G. The experimental works of Kondo et al. [START_REF] Kondo | Nonlinear pressure dependence of the elastic moduli of fused quartz up to 3 GPa[END_REF] have demonstrated that these moduli decrease with increasing pressure up to 2.4 -2.5 GPa, after which they increase with increasing pressure until 3 GPa. Figure 2.15, taken from [START_REF] Kondo | Nonlinear pressure dependence of the elastic moduli of fused quartz up to 3 GPa[END_REF], shows this dependence between the pressure and the elastic moduli. Kondo et al. [START_REF] Kondo | Nonlinear pressure dependence of the elastic moduli of fused quartz up to 3 GPa[END_REF] have also reported that the experimental results associated to this nonlinear behavior are reproducible. It would therefore seem that there are no observable permanent change in these moduli and this behavior is perfectly elastic.

Beyond approximately 3 GPa, the silica glass behavior returns linear elastic until the densification pressure P d ≈ 8 GPa, where the silica glass begins to exhibit signs of permanent deformation. Table 2.1 presents the initial elastic properties of silica glass which are used in the present work.

Permanent deformation: densification

Contrary to crystalline solids, when silica glass, which is an amorphous material, is subjected to high hydrostatic compression (in the range of 8 GPa or beyond), it begins to exhibit signs of permanent deformation which is quite different from the plastic flow observed in crystalline solids. Indeed, the plastic flow in crystalline solids is volume-conservative and initiates under shear stress only. The spherical part of the stress tensor has no influence on the plastic flow. However, the permanent deformation of silica glass is always accompanied by a volume change. In addition, as found in several previous works [START_REF] Bridgman | Effects of Very High Pressures on Glass[END_REF][START_REF] Keryvin | Contribution à l'étude des mécanismes de déformation et de fissuration des verres[END_REF][START_REF] Rouxel | Poisson's Ratio and the Densification of Glass under High Pressure[END_REF], it can initiate under hydrostatic pressure. For these reasons, the permanent deformation of silica glass is generally ). These quantities are linked through the equation of the mass conservation which can be written:

∆V V 0 = - ∆ρ ρ 0 1 + ∆ρ ρ 0 (2.23)
The densification of silica glass which was discovered in the fifties has been widely studied in the literature [START_REF] Brückner | Properties and structure of vitreous silica. I[END_REF][START_REF] Brückner | Properties and structure of vitreous silica[END_REF][START_REF] Bridgman | Effects of Very High Pressures on Glass[END_REF][START_REF] Keryvin | Contribution à l'étude des mécanismes de déformation et de fissuration des verres[END_REF][START_REF] Rouxel | Indentation deformation mechanism in glass: Densification versus shear flow[END_REF]90]. This phenomenon, which initiates at the densification pressure P d , continues to evolve until the saturation pressure P s ≈ 20 GPa and can reach 21 % of density change at this pressure. Beyond the saturation pressure P s , this process stops evolving and the material is reset to its elastic behavior.

As reported in the literature [START_REF] Brückner | Properties and structure of vitreous silica. I[END_REF][START_REF] Brückner | Properties and structure of vitreous silica[END_REF][START_REF] Bridgman | Effects of Very High Pressures on Glass[END_REF][START_REF] Keryvin | Contribution à l'étude des mécanismes de déformation et de fissuration des verres[END_REF][START_REF] Rouxel | Indentation deformation mechanism in glass: Densification versus shear flow[END_REF], during the densification process, the mechanical properties of silica glass show an anomalous dependence on the pressure. Indeed, and contrary to normal glasses, these properties increase with increasing pressure beyond P d . According to Reference [90,[START_REF] Keryvin | Contribution à l'étude des mécanismes de déformation et de fissuration des verres[END_REF], Young's modulus and Poisson's ratio of silica glass present a spectacular augmentation of 46 % and 56 %, respectively, at the end of the densification process (Fig. Figure 2.17: Evolution of silica glass densification with pressure using different samples (ex-situ measures); A: silica powder [START_REF] Cohen | Effects of Ultra high Pressures on Glass[END_REF], B: thin discs [START_REF] Christiansen | Irreversible Compressibility of Silica Glass as a Means of Determining the Distribution of Force in Highpressure Cells[END_REF], C thin discs [START_REF] Bridgman | Effects of Very High Pressures on Glass[END_REF] and D: thick samples [90]. (Taken from [90]) Remark: Several studies measured the densification parameters of silica glass, and different values for each parameter can be found in the literature (Tab. 2.2). Also, these parameters depend hugely on the type of the samples used to study this phenomenon [START_REF] Cohen | Effects of Ultra high Pressures on Glass[END_REF][START_REF] Christiansen | Irreversible Compressibility of Silica Glass as a Means of Determining the Distribution of Force in Highpressure Cells[END_REF][START_REF] Bridgman | Effects of Very High Pressures on Glass[END_REF]90]. For example, the results obtained using powder samples are very different from those obtained using massive samples (Fig. 2.17). The different numerical values given in this paragraph, which have been obtained by Ji [90] using thick samples, are for illustrative purposes only. These values have been used to validate statically the silica glass numerical model developed in the present work.

P d (GPa) P s (GPa) ∆V V 0 (%)
Static studies [START_REF] Brückner | Properties and structure of vitreous silica. I[END_REF][START_REF] Brückner | Properties and structure of vitreous silica[END_REF][START_REF] Keryvin | Contribution à l'étude des mécanismes de déformation et de fissuration des verres[END_REF][START_REF] Bridgman | Effects of Very High Pressures on Glass[END_REF]90] 2 -10 20 -30 7.5 -17.4 Dynamic studies [START_REF] Polian | Sound velocities and refractive index of a-SiO2 to 25 GPa[END_REF][START_REF] Zha | Acoustic velocities and refractive index of SiO2 glass to 57.5 GPa by Brillouin scattering[END_REF][START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF][START_REF] Sugiura | Densified silica glass after shock compression[END_REF][START_REF] Marsh | LASL Shock Hugoniot Data[END_REF][START_REF] Wackerle | Shock wave compression of quartz[END_REF][START_REF] Lalle | Utilisation de l'eau, de la silice fondue et de deux verres de la société Shott comme fenêtre[END_REF] 8.8 -10 16.6 -30 41

Table 2.2: Dispersion range of the silica glass densification parameters

Effects of the shear stresses on the silica glass densification

Almost all the works on the densification of silica glass agree that the densification process can take place under hydrostatic pressure. However the shearing effects on this phenomenon present a central issue and need more investigations. There are several works on this subject in the literature [START_REF] Rouxel | Indentation deformation mechanism in glass: Densification versus shear flow[END_REF][START_REF] Christiansen | Irreversible Compressibility of Silica Glass as a Means of Determining the Distribution of Force in Highpressure Cells[END_REF][START_REF] Mackenzie | High-Pressure Effects on Oxide Glasses: I, Densification in Rigid State[END_REF], some of which are controversial. It began with the works of Bridgman and Šimon [START_REF] Bridgman | Effects of Very High Pressures on Glass[END_REF]. They demonstrated that amorphous silica exhibits permanent densification which translates to a permanent volume decrease of approximately 7.5 % at approximately 20 GPa. They also demonstrated that this propensity decreases when other substances are added, such as sodium oxide, Na2O. However, at that time, technical difficulties prevented the generation of strictly hydrostatic conditions at pressures above 5 GPa [START_REF] Bridgman | Effects of Very High Pressures on Glass[END_REF]. Uniaxial pressures were applied on thin disks to study the densification of glasses. Unfortunately, these tests induced additional shear stresses, which explain the systematic fracturing of the samples. Cohen and Roy [START_REF] Cohen | Effects of Ultra high Pressures on Glass[END_REF] studied the densification of silica powder under uniaxial pressures. They found that the densification begins at approximately 2 GPa (10 GPa in Bridgman and Šimon's work [START_REF] Bridgman | Effects of Very High Pressures on Glass[END_REF]) and the volume reduction can far exceed 7.5 %. According to Reference [START_REF] Weir | Comments on "Effects of Ultrahigh Pressures on Glass[END_REF], the difference between the previous results is due to the shearing effects, which are most important in the case of silica powder and allow a greater densification level to be reached. This explanation was belied later by Cohen and Roy [START_REF] Cohen | Reply to "Comments on'Effects of Ultrahigh Pressures on Glass[END_REF], who excluded the shearing effects on the densification level. For more details, the reader can refer to [START_REF] Keryvin | Contribution à l'étude des mécanismes de déformation et de fissuration des verres[END_REF]. Despite the numerous works on this subject, at present, the role of the shearing stresses in the densification remains ambiguous. Chapter 5 will try to corroborate these works numerically.

Fracture

It is known that the silica glass is the most common example of amorphous brittle materials. The fracture of this material is an especially challenging problem [START_REF] Wiederhorn | Fracture of silicate glasses: ductile or brittle?[END_REF]. It depends on several parameters such as the loading and micro-defects and flaws existing in the material. Also, it is demonstrated that this material can interact with its environment to modify the fracture properties [START_REF] Michalske | A molecular mechanism for stress corrosion in vitreous silica[END_REF][START_REF] Wiederhorn | Influence of Water Vapor on Crack Propagation in Soda-Lime Glass[END_REF]. These properties must be separated between "slow crack growth" and sudden fracture. In the former, the environment, especially water, interacts with strained bonds at the tip of existing cracks and can induce physi-and chemi-sorption processes which weaken the bond and allow the crack to advance slowly. In the latter, the crack is assumed to propagate without interaction with the environment. For the sake of simplicity, the present work only considers the case of sudden fracture in which the failure process is totally determined by mechanical considerations. As will be seen in Section 2.2.2.3, this simplifying assumption might be used without remorse to study the spalling phenomenon which occurs inside the material where the flaws density is very weak. According to Reference [START_REF] Lawn | Fracture of brittle solids[END_REF], brittle fracture is thought to be initiated under tensile stress in mode I. The theoretical strength σ th f of a material without preexisting flaws, that fails by brittle fracture, can be calculated based on a knowledge of its structure and the properties of its inter-atomic bonds. This value is known as the cohesive strength of the material, and it can be calculated for a material given the Young's modulus E, the surface energy γ, the atomic structure and an assumption of the form of the inter-atomic potential [START_REF] Orowan | Fracture and strength of solids[END_REF]. Using the Morse inter-atomic potential function, σ th f can be expressed:

σ th f = E γ 4 r 0 (2.24)
with r 0 is the inter-atomic distance. Applying (2.24) for silica glass, one obtains σ th f = 23 GPa. This value widely exceeds the measured strength σ meas f which is about a few tens of MPa. The Discrepancy between the theoretical and measured strengths is due to the surface flaws. To correct this problem, Griffith [START_REF] Griffith | The Phenomena of Rupture and Flow in Solids[END_REF] has proposed another equation to calculate the fracture strength from a material containing a crack of length c:

σ G f = 2 E γ π c (2.25)
Results obtained using this last Equation (2.25) compare favorably with measured ones, provided that the crack length c is well chosen. However, the choice of this parameter c is not a straightforward issue and represents, at present, an open question. In practice, it is determined experimentally, but the great dispersion of the experimental results makes this way ineffective.

Silica glass response under shock compression

As seen in Section 2.2.1, in statics, the silica glass mechanical behavior can be divided into three domains: elastic domain between 0 and 8 GPa; densification domain between 8 GPa and 20 GPa; saturation domain beyond 20 GPa. Also, it is shown that the elastic domain of this material can be divided into two parts: nonlinear elasticity up to 3 GPa and linear elasticity between 3 and 8 GPa.

In dynamics, these different phenomena have also been reported by several researchers [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF][START_REF] Sugiura | Densified silica glass after shock compression[END_REF][START_REF] Marsh | LASL Shock Hugoniot Data[END_REF][START_REF] Wackerle | Shock wave compression of quartz[END_REF][START_REF] Lalle | Utilisation de l'eau, de la silice fondue et de deux verres de la société Shott comme fenêtre[END_REF]. However, the densification parameters obtained in this regime are different from those obtained in statics (Tab. 2.2). The discrepancies between the static and dynamic results are not well explained in the literature. They can be due to dynamic effects of silica glass or also simply to the dispersion of results as encountered in statics. Especially, analyzing the dynamic studies of silica glass densification, the first point to be noted is the spread of the experimental results which differ from a work to another. For example, the HEL (P d in the case of hydrostatic pressure) varies from 8 GPa [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF] to 10 GPa [START_REF] Lalle | Utilisation de l'eau, de la silice fondue et de deux verres de la société Shott comme fenêtre[END_REF]. Table 2.2 summaries the dispersion of the densification parameters in dynamics. In the present work, the silica glass dynamic behavior will be studied numerically using the densification parameters obtained by Marsh [START_REF] Marsh | LASL Shock Hugoniot Data[END_REF].

Hereafter, the effects of these phenomena on the shock compression response of silica glass will be reviewed.

Propagation of compression wave in the region of reversibility

As a consequence of the nonlinear elastic behavior under approximately 3 GPa, the wave propagation velocity decreases with increasing particle velocity up to approximately 300 ms -1 [START_REF] Wackerle | Shock wave compression of quartz[END_REF][START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF]. This leads to the formation of an elastic precursor which spreads out as the wave propagates into the material. Consequently, in the case of low stresses (under 3 GPa), it is impossible to generate a shock wave in silica glass as reported in several experimental works (Fig. 2.18). Beyond this value, since the silica behavior returns linear elastic, the compression front will steepen up as the wave propagates into the material and a shock wave will be generated. Figure 2.18 shows the particle velocity profiles for different compression levels (taken from [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF]). As can be seen in this figure, there is a sudden jump of particle velocity beyond 3 GPa which indicates the formation of a shock wave.

Propagation of compression wave in the region of irreversibly

Another peculiarity of the silica glass response under shock is the formation of a dual wave in the region of high pressures due to the change in the slope of the Hugoniot curve as mentioned in almost all the previous works studying the silica glass dynamic behavior. The dual wave in the silica glass is very similar to that obtained in the case of elastoplastic solids: the shock front splits into reversible (elastic) compression wave and irreversible (densification) compression wave (Fig.

2.18).

Also, as reported by Lalle [START_REF] Lalle | Utilisation de l'eau, de la silice fondue et de deux verres de la société Shott comme fenêtre[END_REF], the silica glass exhibits a gradual acceleration of the particle velocity beyond the HEL. This can be explained by the augmentation of the bulk modulus K in the densification region. According to the same author [START_REF] Lalle | Utilisation de l'eau, de la silice fondue et de deux verres de la société Shott comme fenêtre[END_REF], this modulus K can reach 400 GPa at the end of the densification process (Fig. 2

Spalling strength

When a non-sustained shock wave propagates into a solid, it is possible to produce high levels of tensile state inside it which can lead to solid spalling. The first point to note is that the spalling fracture strength is generally much greater (in the range of GPa) than the fracture strength measured in quasi-static tests (a few tens of MPa). This can be explained by the fact that the spalling process occurs inside the material where the flaws density is very weak. This process is therefore generated through the internal texture of the material. In contrary, in the case of quasi-static tests, the fracture results from the propagation of the preexisting surface flaws which require lower energy to be activated. In addition, as reported by de Resseguier [START_REF] De Resseguier | Etude expérimentale et numérique du comportement de matériaux de type verre soumis à des ondes de choc[END_REF], the spalling strength depends on the loading velocity, the amplitude and duration of the induced tension and the distance between the tensile zone and the free surface. The smaller this distance, the smaller the spalling strength (activation of the surface defects).

It should be noted that this parameter is very difficult to be experimentally determined. This is attributed to the existence of "failure waves" [START_REF] Bourne | On the shock induced failure of brittle solids[END_REF] immediately behind the shock front and damage the material before the spalling process takes place. Depending on the type of glasses, this parameter variates from 500 MPa for borosilicate glass to more than 2.5 GPa for silica glass [START_REF] Michel | Phénomène d'impact à haute vitesse sur cibles minces fragiles : application au projet de laser mégajoule et à la problématique des débris spatiaux[END_REF].

Summary of the silica glass response under shock compression

The main specificities of the silica glass mechanical behavior and the associated numerical values which are used in the present dissertation can be summarized as follows.

• Under static hydrostatic compression:

-Nonlinear elastic behavior up to 3 GPa due to the pressure dependence of the bulk and shear moduli.

-Linear elastic behavior between 3 GPa and the P d = 8 GPa.

-Densification behavior between the P d and P s = 20 GPa.

-Young's and shear moduli increase respectively by 46 % and 56 % at the end of densification (ex-situ measures).

-Beyond P s , the saturation domain takes place.

-Fracture strength is about a few tens of MPa. In this work, it is fixed at σ f = 64 MPa.

• Under shock compression:

-Nonlinear elastic behavior up to 3 GPa due to the pressure dependence of the bulk and shear moduli.

-Linear elastic behavior between 3 GPa and the HEL = 10 GPa.

-Densification behavior between the HEL and 30 GPa.

-Increase of bulk modulus with densification to 400 GPa at the end of densification process.

-Beyond 30 GPa, the saturation domain takes place.

-Spalling threshold of silica glass is higher than 2.5 GPa. For the qualitative study in Chapter 6, this parameter is fixed at 2.5 GPa.

Numerical simulation

Numerical simulation has increasingly become a very important approach for solving complex practical problems in engineering and science. With the help of the increasing computers power, this modern numerical approach attacks the original problems in all its detail without making too many assumptions, rather than adopting the traditional theoretical practice using assumptions and approximations to simplify the studied problem. It plays a valuable role in providing a validation for theories, offers insights to the experimental results and assists in the interpretation or even the discovery of new phenomena. Numerous numerical methods can be found in the literature which are used to simulate the mechanical behavior of materials, each one has its advantages and drawbacks. These numerical methods can be divided into two main classes, continuum and discrete methods, which follow similar procedure to serve a practical purpose (Fig. 2.22). From the observed physical phenomena, mathematical models are derived with some possible simplifications and assumptions. These mathematical models are generally expressed in terms of governing equations with boundary conditions (BCs) and/or initial conditions (ICs).

To solve numerically the governing equations, the problem domain needs to be discretized, if not naturally discrete, into finite discrete components. Then, the governing equations of the discretized problem and the associated numerical algorithms have to be established. The implementation of a numerical simulation involves translating the domain decomposition and numerical algorithms into a computer code in some programming languages.

Continuum methods

The continuum methods (CMs), such as the Finite Element Method (FEM) [START_REF] Zienkiewicz | Finite Element Method: Its Basis & Fundamentals[END_REF][START_REF] Zienkiewicz | The Finite Element Method for Solid and Structural Mechanics[END_REF][START_REF] Zienkiewicz | The Finite Element Method for Fluid Dynamics[END_REF], the Smoothed Particle Hydrodynamics (SPH) [START_REF] Randles | Smoothed particle hydrodynamics: Some recent improvements and applications[END_REF][START_REF] Liu | Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments[END_REF][START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF], and the Natural Element Method (NEM) [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF][START_REF] Chinesta | Natural element method for the simulation of structures and processes[END_REF][START_REF] Yvonnet | The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces[END_REF][START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous problems: the constrained natural element method (C-NEM)[END_REF][START_REF] Illoul | On some aspects of the CNEM implementation in 3D in order to simulate high speed machining or shearing[END_REF][START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF][START_REF] Sukumar | Natural neighbour Galerkin methods[END_REF], are the most commonly used methods in Computational Mechanics and have been widely applied to various area of Computational Solid Mechanics (CSM) and Computational Fluid Dynamics (CFD). They consist of solving, by means of numerical algorithms, the governing equations stemming from the continuum mechanics which states that matter is a continuum and completely fills the space it occupies. As a first step, the continuum problem domain is approximated by a finite number of discrete components made up of reference points or nodes. Then, the original governing equations are changed into a system of algebraic equations. The resolution of this system gives the solutions at the reference points. For other points in the occupied space, the field variables can be approximated by interpolation or averaging (in the case of SPH) of the solutions at the reference points. The CMs are well adapted for simulating physical phenomena in which the continuity assumption is valid and remains valid during the simulation. They are suited to study phenomena having length scales much greater than the inter-atomic distance. However, this class of method is less suitable to study problems with occurrence of discontinuities in the studied domain, such as wear, fracture and abrasion problems. Indeed, additional treatments must be considered to deal with the new discontinuity surfaces. Depending on the continuum method, there are several techniques used in the literature to circumvent the discontinuities problem.

For the case of the grid-based methods and particularity the FEM, re-meshing or rezoning technique is proposed to take into account the new emerged discontinuity surfaces [START_REF] Gerstle | Finite and boundary element modeling of crack propagation in two-and three dimensions[END_REF][START_REF] Gerstle | Three-Dimensional fatigue crack propagation analysis using the boundary element method[END_REF][START_REF] Martha | Arbitrary crack representation using solid modeling[END_REF]. This technique, which involves overlaying of a new mesh on the old one, can be tedious and very time-consuming. Moreover, since the field variables in the new mesh nodes are approximated from the old ones, some approximation errors can be introduced and the numerical simulation loses accuracy with each re-meshing. Another technique widely applied in linear elastic fracture analysis [START_REF] Griffith | The Phenomena of Rupture and Flow in Solids[END_REF] is to use singular elements [START_REF] Henshell | Crack tip nite elements are unnecessary[END_REF][START_REF] Barsoum | Triangular quarter-point elements as elastic and perfectlyplastic crack tip elements[END_REF]2].

In spite of the successes using finite elements in computational fracture, mesh generation in three dimensions is time-consuming and especially burdensome for multiple crack configurations and crack growth simulations. It is difficult to explicitly model the crack topology as part of the finite element since accuracy considerations require significant refinement in the vicinity of the crack front. More recently, an extension of FEM, called the Extended Finite Element Method (X-FEM), is developed by [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] to alleviate the FEM shortcomings associated with meshing of discontinuities surfaces. In this approach, the standard displacement-based finite element approximation is enriched by additional (special) functions using the framework of unity partition. This method is the most suitable grid-based method to model evolving discontinuities such as crack growth. However, the use of this method to model problems with muti-cracks remains problematic. The difficulties and limitations of this method and the grid-based methods in general are especially evident when simulating hydrodynamic phenomena such as High Velocity Impact (HVI) where multi-cracks can be generated.

A recent strong interest is focused on the development of the next generation of continuum methods, meshfree methods, such as SPH [START_REF] Randles | Smoothed particle hydrodynamics: Some recent improvements and applications[END_REF][START_REF] Liu | Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments[END_REF][START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF], Finite Point Method (FPM) [START_REF] Oñate | A Stabilized Finite Point Method for Analysis of Fluid Mechanics Problems[END_REF][START_REF] Oñate | A Mesh-Free Finite Point Method for Advective-Diffusive Transport and Fluid Flow Problems[END_REF][START_REF] Oñate | A Finite Point Method for Incompressible Flow Problems[END_REF] and Element Free Galarkin Method (EFGM) [START_REF] Belytschko | Element-Free Galerkin methods[END_REF]. The key idea of these methods is to provide accurate and stable numerical solutions for the governing equations with a set of arbitrarily distributed nodes (or particles) without using any mesh or connectivity between them. Since the computational frame in the meshfree methods is a set of arbitrarily distributed nodes rather than a system of pre-defined mesh/grid, these methods are attractive in dealing with problems that are difficult for traditional grid-based methods. This is despite the fact that they are very time-consuming and suffer from some numerical problems such as boundary inaccuracy and difficulties to impose essential (Dirichlet) boundary conditions. In addition, they produce results with errors larger than those obtained using grid-based methods.

Discrete methods

The discrete methods (DMs), such as the Distinct (or Discrete) Element Method (DEM) [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], the Generalized Discrete Element Method (GDEM) [START_REF] Williams | The Theoretical Basis of the Discrete Element Method[END_REF], and the Discontinuous Deformation Analysis (DDA) [START_REF] Shi | Discontinous deformation Analysis[END_REF], are a family of numerical methods consisting in simulation of a set of rigid (or pseudo-deformable) bodies in interaction with each other by contact laws. The Fundamental Principle of Dynamics (FPD) is applied to compute the solids displacements and rotations. Originally, this class of methods has been developed for problems of granular materials in rocks mechanics [START_REF] Cundall | Computer Model for Simulating Progressive Large Scale Movements in Blocky Rock Systems[END_REF]. More recently, this class of methods has been applied to study the damage of heterogeneous media such as concrete [START_REF] Hentz | Discrete element modelling of concrete submitted to dynamic loading at high strain rates[END_REF] and rocks [START_REF] Bobet | Numerical models in discontinuous media: review of advances for rock mechanics applications[END_REF] , but also to study the damage of homogeneous media such as ceramics [START_REF] Tan | Discrete element method (DEM) modelling of fracture and damage in the machining process of polycrystalline sic[END_REF] and glasses [4,3]. In these studies, the material is modeled by an agglomerate of discrete elements which interact via bilateral cohesive links to ensure the cohesion of the medium. The cohesive links (Fig. • choice of the cohesive links and identification of their microscopic parameters to ensure a correct macroscopic behavior

• construction of DEM domains taking into account the structural properties of the problem domain, e.g. homogeneity and isotropy, such that the macroscopic mechanical properties are independent of the discrete elements number.

Several works can be found in the literature which propose solutions of these points. Among them, one can cite the recent works of André et al. [4]. Nowadays, DMs present an alternative way to study physical phenomena requiring a very small scale analyses, where the continuity assumption is no longer valid, or problems with discontinuities that cannot be easily treated by CMs. Indeed, in DMs, a continuum is divided by fully persistent discontinuities delimiting a finite number of particles interconnected through these discontinuities. Consequently, this class of methods has become increasingly used to study the cracking behavior of continuum materials such as silica glass [4,3]. However, DMs are very costly in terms of CPU which limits the length of the DM simulation and/or the number of particles used to model the continuum. Therefore, they are only applied to model small domains.

Coupling methods

As seen above, there are several numerical methods used in Computational Mechanics which can be classified into continuum and discrete methods according to the physics behind them. Each numerical method has its features and advantages which can be complementary. This suggests that it would be computationally beneficial to combine two or more methods so as to strengthen their advantages and to avoid their drawbacks. This idea has led to the development of several coupling methods between continuum methods (continuum-continuum coupling) or between continuum and discrete methods (discrete-continuum coupling). These coupling approaches consist in applying different numerical methods together in separate regions (or with overlap) in the problem domain (Fig.

2.25).

This technique is widely used in the literature to study multi-scale phenomena, in which different methods are recommended for each scale, and to "zoom in" on some regions in the problem domain where complex phenomena (e.g. cracking) can take place. The continuum-continuum coupling is generally applied to couple a mesh-free method and a grid-based one. The former, which is generally costly but treats better the discontinuity problems, is applied where cracks can appear; the later is applied in the rest of the problem domain. Concerning the discrete-continuum coupling, the discrete method which is very costly is applied only in the regions in which complex phenomena can take place, and the continuum method is applied in the rest of the problem domain to reduce the computation time.

What class of numerical methods best meets the objectives of this dissertation?

The main objective of this dissertation is to simulate the laser-induced damage in silica glass. This highly dynamic application includes several complex phenomena at different scales. Moreover, discontinuities problems can occur in some regions due to cracking. Therefore, the continuum grid-based methods, which cannot properly treat discontinuities, will be excluded from the numerical methods to be used in the loading regions. Also, the continuum mesh-free methods, which suffer from several difficulties related to the boundaries treatment, are not a good choice to model these regions. Consequently, even the discrete methods are very time-consuming, they represent an obvious candidate to model the complex phenomena in the loading regions. However, in most situations and particularly in laser applications, the effects requiring very small scale analyses by discrete methods are localized in small portions of the studied domain. Thus, the use of a discrete-continuum coupling method to treat the phenomena at each scale seems to be advantageous. The discrete method will be applied only in the high loaded regions, whereas the continuum method will be applied in the rest of the studied domain. Therefore, a new discrete-continuum coupling approach between a discrete method and a continuum method will be developed in this dissertation. The choice of the two numerical methods to be coupled will be discussed later in Chapter 3.

Conclusion

This chapter has reviewed the principle scientific advancements in the different disciplines interacting to accomplish the goals of this dissertation which is the study of the laser-induced damage in silica glass. In Section 2.1, a basic knowledge of the laser-matter interaction theory is given. The physics behind this theory can be divided into two steps: generation of shock wave by laser ablation and shock wave propagation within the irradiated material. The first step is not within the scope of this work which aims to study the mechanical behavior of silica glass under laser shock processing. The temporal pressure profile generated by the laser-matter interaction is assumed to be given and will directly be applied on the shock surface (irradiated surface). During its propagation within a material, the shock wave is subject to several mechanical phenomena such as hydrodynamic damping, elastoplastic damping or also reflexion and transmission (when migrating between media). In addition, it can damage the considered material. These different phenomena will be highlighted through the qualitative mechanical study of the Laser Shock Processing on silica glass.

In Section 2.2, the complex mechanical behavior of silica glass in described. This material which belongs to the category of anomalous glasses shows an anomalous behavior in its thermal and mechanical properties. In the region of ordinary pressures, it behaves in a perfectly elastic manner up to the densification pressure P d . Beyond this value, it begins to exhibit signs of permanent deformation called densification up to the saturation pressure (P s ) where this material returns elastic. In addition, it presents a quite different cracking behavior compared to other glasses. All these specificities will be taken into account to model the silica glass mechanical behavior in Chapter 5.

Finally, a brief description of the different numerical methods used to simulate mechanical problems is given in Section 2.3. These methods can be divided into two classes: continuum methods and discrete methods. The first class is widely used in Computational Mechanics to study coarse scale phenomena having length scales much greater than the inter-atomic distance. However, it presents several difficulties to study complex phenomena requiring very fine scale analyses or phenomena with discontinuities. The second class presents an alternative way to study these complex phenomena. But, it is generally very time-consuming which limits its application to small domains compared to the analysis scale. Since the LSP applications require a very small scale analysis in the shocked regions to capture the microscopic effects due to the shock wave propagation, it is recommended to use a discrete method in these regions. Whereas, a continuum method can be used outside of these regions to reduce the computational cost. Therefore, a discrete-continuum coupling approach adapted for dynamic simulations is proposed in this dissertation to simulate the LSP processing on silica glass. The discrete and continuum methods involved in this approach will be chosen in the next chapter.

Introduction

In the literature, there are several numerical methods used in Computational Mechanics to simulate the mechanical behavior of materials. As shown in the previous Chapter, these methods can be classified into two main classes: discrete methods and continuum methods. The first class is based on the Newtonian mechanics and attempts to model a material by an agglomerate of discrete elements. This facilitates the treatment of discontinuities which may occur during simulations, but it is very time-consuming. The second class is based on the continuum mechanics which considers a material as a continuum. Simulations using this class of methods are generally faster and less expensive compared to the first class, however, additional treatments need to be performed to deal with complex phenomena involving discontinuities such as fracture. The features and advantages of these two classes are largely complementary. Therefore, it would be advantageous to combine two methods from these classes to strengthen their advantages and circumvent their drawbacks. The question that arises here is how to choose the numerical methods which will be coupled in this work to simulate the Laser Shock Processing (LSP) on silica glass. This chapter aims to answer this question.

As can be seen from the literature, each class includes numerous numerical methods, some of which share several fundamental concepts. Therefore, to simplify the choice of the appropriate numerical methods, it would be helpful to divide each class into several subclasses, each one including the methods that share the main features in their principles. A first classification involving only subclasses (not methods) will be made, according to their advantages and drawbacks. This will allow us to retain the discrete and continuum subclasses that best comply with the objectives of this work. The appropriate methods will then be selected from the retained subclasses only.

Choice of the discrete method

The different discrete methods found in the literature can be divided into three main categories (subclasses) as shown in Figure 3.1. The fundamental concepts of each one will be briefly recalled hereafter.

Lattice models

In Lattice models, a solid is modeled by a set of nodes connected with truss or beam elements [START_REF] Schlangen | Experimental and numerical analysis of micromechanisms of fracture of cement-based composites[END_REF][START_REF] Schlangen | Simple lattice model for numerical simulation of fracture of concrete materials and structures[END_REF]. Typically, nodes have neither masses nor volumes (they do not occupy volumes). The resolution of mechanical problems with this class of discrete methods is based on the construction of a global stiffness matrix from the local Figure 3.1: Classification of discrete methods connection properties. Both regular and irregular lattices were studied. Originally, the lattice was used to represent elastic continuum; the equivalence was established for both truss [START_REF] Hrennikoff | Solution of problems of elasticity by the frame-work method[END_REF] and beam [START_REF] Schlangen | New method for simulating fracture using an elastically uniform random geometry lattice[END_REF] elements. Later on, obvious enhancements, such as brittle beam failure, were introduced. Lattice models nicely show the emergence of relatively complex structural behaviors, although fairly simple formulas are used to describe the governing local processes.

The major drawback of these models is that the nodes do not have volumes, which can cause numerical problems related to crack closure in post-fracture stage. Solutions are given to circumvent this problem. These solutions consist of associating to each element an equivalent volume, based on the spacial Voronoï decomposition [START_REF] Ragueneau | Inelastic behavior modelling of concrete in low and high strain rate dynamics[END_REF] for example.

Particle models

This class of methods is very close to the first discrete approach proposed in the literature by Cundall and Strack [START_REF] Cundall | Computer Model for Simulating Progressive Large Scale Movements in Blocky Rock Systems[END_REF][START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]: the discrete (distinct) element method (DEM). Contrary to Lattice models, particle models consider elements with masses and volumes. These elements often have disc shape (in 2D) or spherical shape (in 3D): only a single parameter (the radius) is required to determine the elements geometry and there is only one possible contact easily detectable between them. Consequently, computer memory requirements and processing time are minimized with these element shapes, even when a relatively large number of elements is used. Nevertheless, discs and spheres can roll or rotate easily. This does not reflect the expected behavior for several materials, in the case of large shear processes for example. To solve this problem, more complex shapes such as ellipses [START_REF] Ting | An ellipse-based discrete element model for granular materials[END_REF], ellipsoids [START_REF] Lin | A three-dimensional discrete element model using arrays of ellipsoids[END_REF], polygons [START_REF] Issa | Numerical analysis of micromechanical behaviour of granular materials[END_REF] and polyhedra [START_REF] Cundall | Formulation of a three-dimensional distinct element model -Part I, A scheme to detect and represent contacts in a system composed of many polyhedral blocks[END_REF] were proposed in the literature to provide more flexibility in element characterization in particle models.

Basically, the associated algorithm involves two stages. In the first stage, interaction forces are computed when elements slightly interpenetrate each other. This force-interpenetration formulation is generally referred to as a "Smooth contact" method or "Force-displacement" method. Actually, the interpenetration between discrete elements, which makes no mechanical sense, represents the relative deformation of the elements surface layers. In the second stage, Newton's second law is applied to determine the acceleration of each element, which is then integrated to find the new element velocity and position. This process is repeated until the simulation is achieved. The associated governing equations are solved using "dynamic explicit" schemes.

Contact dynamics

This class of methods provides an alternative approach based on a "non-smooth" formulation of mutual exclusion and dry friction between elements [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF][START_REF] Luding | Simulations of pattern formation in vibrated granular media[END_REF]. It introduces the notion of non-smooth (irregular) contact between elements which is, at present, subject of several studies. Interpenetration between elements is prevented: no elastic contact law is used between them. There are two main numerical integrators for these methods, which are of "dynamic implicit" type: the event-driven integrators, also referred to as the Even-Driven Method (EDM) [START_REF] Luding | Simulations of pattern formation in vibrated granular media[END_REF], and the socalled time-stepping integrators, also referred to as the Contact Dynamics Method (CDM) [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF].

In EDM, a collision or "event" occurs when two rigid elements touch each other and the post-collisional and angular velocities are prescribed by a collision operator [START_REF] Rapaport | The event scheduling problem in molecular dynamic simulation[END_REF]. Despite being very accurate, the even-driven integrators treat only one force at a time. Therefore, they are not well adapted for problems with many simultaneous contacts, as often encountered in mechanics.

To overcome this limitation, Jean and Moreau [START_REF] Jean | The non smooth contact dynamics method[END_REF][START_REF] Moreau | Nonsmooth mechanics and applications[END_REF] have developed the Contact Dynamic Method (CDM) which has a specialized numerical scheme for problems with many contacts. The governing equations are expressed as differential inclusions and the accelerations are replaced by velocity jumps. In the generic CDM algorithm, an iterative process is used to compute forces and velocities. It consists of solving a single contact problem with all other contacts kept constant, and iteratively updating the forces until a convergence criterion is fulfilled. Two basic kinematic constraints are used between elements in the CDM formulation:

• The Signorini conditions which state that the normal force f n is repulsive when the elements are in contact (distance between them is zero), and f n = 0 otherwise. Besides, a kinematic constraint is added for persistent contact: f n vanishes if the relative velocity between elements in contact is zero.

• The Coulomb's friction law, which relates the sliding velocity and the friction force f t .

These kinematic constraints can also be completed with a "rolling friction" which introduces a moment resistance [START_REF] Bratberg | Dynamic rearrangements and packing regimes in randomly deposited two-dimensional granular beds[END_REF].

Within the scope of CDM, the time resolution is much larger than the collision characteristic time (unlike the particle approaches). Instead, the time step represents a unit of time during which collisions can occur causing velocity jumps. Although CDM has been successfully used for several geomechanical problems [START_REF] Donzé | Advances in Discrete Element Method applied to Soil, Rock and Concrete Mechanics[END_REF], it is much more difficult to implement than the Lattice and particle models. Also, the prediction of the contact forces and particle velocities in the following step from the current configuration presents a challenge and is currently the subject of several studies.

Classification and choice of the discrete method: DEM

The Contact Dynamics models are generally used to study quasi-static problems or problems with relatively low dynamic effects. This class of methods is perfectly suitable to study granular materials. However, the LSP processing on silica glass is a highly dynamic problem. Besides, the silica glass is a brittle elastic material. Therefore, the use of models based on regular or "smooth" interaction laws seems to be advantageous since the elasticity is naturally taken into account by these interaction laws. The Contact Dynamics models will thus be disqualified from the competition.

The features of the two other subclasses (Lattice and particle models) are largely complementary. Indeed, in Lattice models, additional treatments must be made to deal with fracture problems (crack closure). This problem is not encountered in particle models in which elements have proper volumes. However, the latter cannot correctly model a continuum using simply disc or spherical elements, especially when shear effects are significant. This problem can be solved using cohesive beams between elements, such as in Lattice models. The complementary advantages of these models has given rise to a family of "hybrid" methods which are between Particle and Lattice models, i.e. by considering sphere elements connected by cohesive beams [3,4,[START_REF] Griffiths | Modelling of elastic continua using a grillage of structural elements based on discrete element concepts[END_REF].

The choice of an "hybrid" model, allowing to efficiently combine the advantages of Lattice and Particle models, is thereby justified. Specifically, the variant of the Discrete Element Method (DEM) recently developed by André et al. [3,4] is retained in this work. This method allows to model an elastic solid by a set of discrete elements linked by cohesive beams; the equivalence has been demonstrated and well explained by André et al. [3,4]. The main specificities and features of this method will be detailed later in Chapter 4.

Choice of the continuum method

As can be noted in the literature, numerous continuum methods are used to simulate the mechanical behavior of materials. Each method has its advantages and drawbacks. These methods can mainly be classified into two subclasses: grid-based methods and meshless methods (Fig. 3.2). This section will review the main features of these subclasses. As a result of this review, an appropriate continuum method will be retained to be coupled with the Discrete Element Method (DEM).

Grid-based methods

In grid-based methods the spacial domain is often represented by discretized elements or meshes. These elements are called "meshes" in the case of Finite Element Method (FEM) [START_REF] Zienkiewicz | Finite Element Method: Its Basis & Fundamentals[END_REF], grids in the case of Finite Difference Method (FDM) [START_REF] Smith | Numerical Solution of Partial Differential Equations: Finite Difference Methods[END_REF] and volumes or cells in the case of Finite Volume Method (FVM) [START_REF] Versteeg | An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Pearson education limited[END_REF]. The terminologies of grids, volumes, cells, and elements carry different physical meanings linked to physical problems. However, all of them can be termed meshes according to the following definition. A mesh, in the broadest sense of the word, is defined as an open space or interstice between the strands of a network which is formed by connecting nodes in a predefined manner. The key here is that, in grid-based methods, mesh must be predefined to provide a certain relationship between the nodes. This Based on a properly predefined mesh, complex ordinary differential or partial differential governing equations can be approximated by a set of algebraic equations for each discretized element (mesh). The system of algebraic equations for the whole problem domain can be formed by assembling the elementary algebraic equations for all the meshes.

There are two main frames for describing the governing equations of mechanical problems: Lagrangian description and Eulerian description. The former, also called material description, is a way to look at material motion where the observer follows an individual material point as it moves through the space and time. The latter, also called spatial description, is another way to look at material motion where the observer focuses on a specific space position through which the material passes over time. The use of one description or the other leads to quite different governing equations. The difference between them arises from the definition of the total time derivative as the local time derivative and the convective derivative in the case of Eulerian description (3.1):

d dt total time derivative = ∂ ∂t local time derivative + V .∇ convective derivative (3.1)
with V represents the velocity and ∇ denotes the vector differential operator (gradient).

Lagrangian methods

Lagrangian methods which are typically represented by the FEM [START_REF] Gouri | Une présentation de la méthode des éléments finis[END_REF][START_REF] Zienkiewicz | Finite Element Method: Its Basis & Fundamentals[END_REF] are based on the Lagrangian description. In this class of methods, the grid or mesh is attached to the material in the entire computation process, and moves with it (Fig. 3.

3).

Each node in the mesh follows the material path at the grid point. As a result, the mesh elements can be deformed during the simulation process. When the problem domain deforms, the mesh elements deform accordingly. Mass, momentum and energy are transported with the movement of the mesh. Because the mass within each cell (or element), remains constant, no mass flux is allowed through the mesh cells boundaries. The Lagrangian methods present several advantages:

• It is very easy to treat the boundary conditions at free surfaces, moving boundaries and material interfaces. Indeed, they are automatically imposed, tracked and determined simply by placing some grid nodes along the boundaries and the material interfaces.

• Complex geometries can be conveniently studied using an irregular mesh.

• Meshing is only required within the problem domain (no additional grids are required beyond it, as will be seen for Eulerian methods). Therefore, Lagrangian methods are computationally efficient.

• Since the mesh is attached to the material, the entire time history of all field variables at a specific material point can easily be obtained during the simulation.

• In the Lagrangian description, the total time derivative is the same that the local time derivative. Therefore, no convective terms exist in the associated governing equations. This makes Lagrangian simulations simpler and faster as no computational efforts are required to treat the convective terms.

Due to these advantages, Lagrangian methods have become very popular to simulate Computational Solid Mechanics (CSM) problems, where deformations are relatively small. However, applying these method in solving problems with large deformations presents a huge challenge. Indeed, in such problems, the mesh can become extremely distorted and stability as well as convergence difficulties can arise leading to computation breakdown. The accuracy of the formulation and consequently the solution can be severely affected. Also, the time step, controlled by the size of the smallest element grid, can become too small to be efficient for computation.

A solution to overcome this inconvenience is to re-mesh (rezone) the problem domain or simply the regions where the initial mesh is hugely distorted. This technique involves overlying of a new (undistorted) mesh on the old one. The computation is then resumed on the new mesh. The field variables are approximated at the new grid nodes via mass, momentum and energy transport using an Eulerian description. Despite the popularity of the rezoning techniques to simulate large deformation problems, they suffer from several difficulties. Indeed, rezoning procedure can be tedious and time-consuming. Besides, the transport of the field variables from old to new mesh is generally accompanied by material diffusion which can lead to a loss of the material history [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF].

Eulerian methods

In Eulerian methods, contrary to Lagrangian methods, the mesh is fixed on the space occupied by the simulated object which, for its part, can move across the fixed mesh cells (grids) as shown in Figure 3.4. All the mesh cells as well as the associated nodes remain fixed in the space and do not move with the material moving across the mesh. The flux of mass, momentum and energy across mesh cell boundaries are simulated to compute the distribution of different field variables in the problem domain. As a result, the shape and the volumes of the grid cells remain constant during the entire computation process. Therefore, simulation of problems with large deformations is possible with Eulerian methods. Since, they do not induce any cells deformation and thus do not cause the same kind of numerical problems as in Lagrangian methods. For these reasons, Eulerian methods are dominant in Computational Fluid Dynamics (CFD) problems and problems with large deformations such as explosion, high velocity impacts (HVI), etc.

Despite the great success of Eulerian methods, there are many disadvantages which limit their application, from which one can cite:

• Since the Eulerian methods simulate the flux of mass, momentum and energy across mesh cells boundaries, it is very difficult to get out the time history of field variables at a fixed material point.

• Eulerian methods require meshing not only within the problem domain, but also beyond. The mesh must be large enough to cover the area through which the simulated object can move. Therefore, coarse grid is sometimes recommended for computational efficiency which can affect the solution accuracy.

• The free surfaces, moving and deformable boundaries and moving material interfaces are difficult to be determined accurately.

• It is difficult to treat irregular and complicated geometries of materials. Generally, tedious mesh generation procedure to convert the irregular geometry of problem domain into a regular computational domain is necessary. This can become sufficiently expensive to be efficient. Efforts in this direction are still ongoing.

Combined Lagrangian-Eulerian methods

The features of the Lagrangian and Eulerian methods are largely complementary. This has led to the development of several approaches involving in the same time the Lagrangian and Eulerian descriptions. The Coupled Eulerian Lagrangian (CEL) [START_REF] Xie | Application of a coupled Eulerian-Lagrangian method to simulate interactions between deformable composite structures and compressible multiphase flow[END_REF][START_REF] Mair | Hydrocode methodologies for underwater explosion structure medium/interaction[END_REF][START_REF] Felici | A coupled Eulerian/Lagrangian method for the solution of threedimensinal vortical flows[END_REF] and the Arbitrary Lagrangian Eulerian (ALE) [START_REF] Stein | Encyclopedia of Computational Mechanics[END_REF][START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF] can particularly be cited.

The CEL approach consists in employing both the descriptions in separate (or with overlap) regions of the problem domain. The Lagrangian region continuously interacts with the Eulerian one through a coupling module allowing for exchange of computational information.

The ALE approach is conceptually close to rezoning techniques for Lagrangian methods. In this case, the mesh moves independently of the material so as to minimize the mesh distortion. In the ALE approach, Lagrangian motion is firstly computed at each time step, followed by possible rezoning step in which the mesh is either rezoned to the original shape or rezoned to some more advantageous shape (between the previous shape and the current shape).

The combined Lagrangian-Eulerian methods are very promising and have received much research interest. However, they are generally difficult to implement. Besides, unexpected termination of the computation process is often encountered, due to a highly distorted Lagrangian element or a very small Eulerian cell.

Meshless methods

Despite the great success of the grid-based methods in both CSM and CFD, their disadvantages have pushed researchers to seek for new numerical methods. One important goal of the initial research is to modify the internal structure of the gridbased methods to become more adaptive and more robust. More effort was given to problems which cannot easily be treated by this class of methods, such as problems with moving boundaries (for Eulerian methods) and large deformation (for Lagrangian methods). Such an effort has given rise to the next generation of computational methods: the meshless methods. The main key of the meshless methods is to provide accurate and stable solutions for integral equations or Partial Differential Equations (PDE) with a set of arbitrary distributed nodes (or particles) without any connectivity between them. Contrary to grid-based methods, this class of methods normally do not require any predefined mesh to establish the system of algebraic equations for the whole problem domain. Instead, only a set of nodes scattered within the problem domain and along the boundaries are used to represent the simulated geometry (Fig. 3.5). In such a representation, no connectivity between nodes is required. The ideal requirements of the meshless methods are: • non-necessity of a predefined mesh, at least in field variable interpolation;

• non-necessity of a predefined mesh at all throughout the computation process.

However, the meshless methods developed so far are not really ideal considering these two points. For example, some meshless methods require background cells (mesh) to derive the system of algebraic equations. The Element Free Galerkin (EFG) method [START_REF] Belytschko | Element-Free Galerkin methods[END_REF] and the Meshless Local Petrov-Galerkin (MLPG) method [START_REF] Atluri | A new meshless local Petrov-Galerkin (MPLG) approach in computational mechanics[END_REF], which use background cells for integration of system matrices driven from the weak formulation, belong to this category.

In the literature, one can find a number of meshless methods, some of which share several common features. It should be noted that almost all the meshless methods are based on the Lagrangian description. Therefore, construction of shape functions is necessary before or during the computation process. These shape functions will be next used to approximate the field variables at any point x within the problem domain as follows (3.2):

u(x) = n ∑ i=0 φ i (x)u i (3.2)
where u is the field variable to be approximated (such as displacement), n is the number of all the nodes, u i is the nodal field variable at the node i, and φ i is the shape function associated to the node i. According to the type of shape functions, the meshless methods can be classified into two main categories: approximation methods and interpolation methods. The main features of each category will be detailed hereafter.

Approximation methods

Contrary to Lagrangian grid-based methods, the shape functions in meshless methods are constructed for the nodes (not for element grids). In the case of approximation methods, these functions generally do not satisfy the Kronecker delta function property at the nodes which is defined as (3.3):

φ i (x j ) = 1 i = j 0 i = j (3.3)
In other words, an approximated field variable using these functions, called approximation shape functions, does not necessary pass through the nodal values (u

(x i ) = u i ).
The first method developed in this category is the famous Smoothed Particle Hydrodynamics (SPH) method [START_REF] Randles | Smoothed particle hydrodynamics: Some recent improvements and applications[END_REF][START_REF] Liu | Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments[END_REF][START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF][START_REF] Lucy | Numerical approach to testing the fission hypothesis[END_REF]. This method was originally invented for modeling astrophysical phenomena, and later it is widely extended for application to problems of continuum solid and fluid mechanics [START_REF] Lucy | Numerical approach to testing the fission hypothesis[END_REF]. The SPH method and its different variants, e.g. the Corrective Smoothed Particle Method (CSPM) [START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF] and Discontinuous Smoothed Particle Hydrodynamics (DSPH) [START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF], are the major type of approximation methods, and have been incorporated into many commercial codes. Several other methods belonging to this category can be found in the literature. Among them, one can cite the Generalized Finite Difference Method (GFDM) [START_REF] Liszka | The finite difference method at arbitrary irregular grids and its applications in applied mechanics[END_REF], the Element Free Galerkin Method (EFGM) [START_REF] Belytschko | Element-Free Galerkin methods[END_REF], the Reproducing Kernel Particle Method (RKPM) [START_REF] Li | Meshfree and particle methods and their applications[END_REF], etc. Several advantages have made the approximation methods very attractive:

• Since no predefined grid is required, approximation methods are able to deal with problems involving extremely large deformations.

• Approximation methods are more adapted than the Eulerian grid-based methods to study problems with free surface, deformable boundary and moving interfaces.

• Approximation methods are more adaptive to treat complex geometries which are represented only by a set of scattered nodes with no connectivity between them.

• Applications of this class of methods are very wide, ranging from micro-scale to macro-scale and even to astronomical scale, and from discrete systems to continuum systems.

• Due to intensive research efforts, several approximation methods such as the SPH are quickly approaching their mature stage.

• The approximation methods (meshless methods) are more adapted to study the fracture of materials. Indeed, they do not use any kind of mesh which can affect the cracking mechanisms.

• The accuracy of solutions is very weakly affected by the relative positions of the nodes. Hence, it is very easy to add or remove nodes. This point is particularly attractive to define an adaptive mesh refinement strategy for problems with localization.

However, the approximation methods have several drawbacks related to the imposition of boundary conditions and implementation requirements:

• To impose essential boundary conditions in a straightforward way as in FEM, it is necessary that: (i) the approximated field variables pass through the nodal value (interpolation shape functions) and (ii) the influence of the internal nodes vanishes at the border of the problem domain (linear shape functions at the border). Since, for approximation methods, generally any of these conditions is verified, it is not easy to impose essential boundary conditions. Therefore, special techniques are required to overcome this problem. In the literature, one can find several works on this subject and several techniques have been proposed. But, in most cases, these techniques are tedious and involve additional computational efforts which can limit the application of these methods to problems with complex/moving interfaces and boundaries. Among these techniques, one can cite the Lagrange multiplier method [START_REF] Belytschko | Element-Free Galerkin methods[END_REF], Penalty method [START_REF] Zhu | A modified collocation method and a penalty formulation for enforcing the essential boundary conditions in the element free galerkin methods[END_REF] and the introduction of singular weight functions [START_REF] Kaljevic | An improved element free Galerkin formulation[END_REF].

• The second drawback of the approximation methods is related to the numerical integration. Actually, almost all the shape functions used in meshless methods are rational (not polynomial). This makes the Gauss Integration schemes not well adapted. Dolbow et al. [START_REF] Dolbow | Numerical integration of the Galerkin weak form in meshfree methods[END_REF] have demonstrated that, using these techniques, a great number of Gauss points is required which can induce an extensive additional cost of computation. Concerning the other popular integration technique, Direct Nodal Integration, Chen et al. [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] have demonstrated that the application of this technique involves numerical instability problems.

• Another drawback of this class of methods is related to the definition of the influence area of the particles (nodes), called support domain or also influence domain. In most cases, the support domain of a particle is taken as a sphere or parallelepiped centered on the particle. As discussed in several previous works [START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF][START_REF] Liu | Smoothed Particle Hydrodynamics (SPH): an Overview and Recent Developments[END_REF], this support should contain a sufficient number of neighbor particles for computation stability and then be sufficiently wide. However, a very wide support domain can affect the accuracy of the solutions. The choice of the size of the support domain is not trivial and represents, at present, an open question [START_REF] Lee | On error estimation and adaptive refinement for element free galerkin method : Part I : stress recovery and a posteriori error estimation[END_REF].

• As reported in several previous works [START_REF] Belytschko | Meshless methods : an overview and recent developments[END_REF][START_REF] Lee | On error estimation and adaptive refinement for element free galerkin method : Part I : stress recovery and a posteriori error estimation[END_REF], the application of approximation methods to non-convex geometries, e.g. the presence of concave boundaries, cavity or cracks, is very challenging. This is mainly due to difficulties to define the support domain in the vicinity of the concave borders. In this case, the support domain can contain nodes which are not really neighbors, e.g. nodes on opposite sides of a crack.

• Finally, compared to Lagrangian grid-based methods, approximation methods generally leads to poor and inaccurate solutions. The accuracy problem is much more amplified in the vicinity of the border of the problem domain.

Despite all the advantages of the approximation methods (meshless methods), the different drawbacks of these methods have limited their application to several mechanical problems.

Interpolation methods

The meshless methods were initially developed to deal with complex problems in both CSM and CFD which cannot be easily treated by grid-based methods, such as fracture and large deformation problems. However, almost all the first generation of meshless methods belong to the category of approximation methods. Consequently, several difficulties related to implementation and boundary condition treatments limit their application. This pushed researchers to make more effort to overcome these difficulties. Such an effort has led to the development of a next generation of meshless methods: interpolation methods. In this generation, the field variables are approximated using shape functions verifying the Kronecker delta function property at the nodes (3.3) (interpolation shape functions). In other words, the approximated field variables pass through the nodal values (u

(x i ) = u i ).
The first method developed in this category is the Natural Element Method (NEM) [START_REF] Braun | A numerical method for solving partial differential equations on highly irregular evolving grids[END_REF][START_REF] Traversoni | Natural neighbour finite elements[END_REF]. In this method, the interpolation shape functions are constructed based on the notion of Natural Neighbor (NN) interpolation. The NN interpolation is defined as a multivariate data interpolation scheme [START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF], which has primarily been used in data interpolation and modeling of geophysical phenomena. This relies on the concepts of Voronoï diagrams [START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques[END_REF] and Delaunay tessellations [START_REF] Delaunay | Sur la sphère vide. A la mémoire de Georges Voronoï[END_REF], which are widely used in computational geometry, to construct the shape functions (called NN shape functions). This method has overcome most difficulties encountered with the earlier generation of meshless methods, such as imposition of Dirichlet boundary conditions and definition of the support domain of the nodes. However, the NN shape functions used in NEM are not strictly linear over non-convex boundaries which is essential for direct imposition of Dirichlet boundary conditions (the test functions do not vanish over the whole region of essential boundaries). This is because NEM is closely related to the convex hull of the given set of points. Roughly speaking, a solution may be found by generalizing the concept of convex hull to non-convex domains.

Recently, Cueto et al. [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF][START_REF] Cueto | Overview and recent advances in Natural Neighbour Galerkin Methods[END_REF] have proposed an extension of NEM, called α-NEM, based on the concept of α-shapes or α-complexes, which are widely used in the field of scientific visualization. These techniques allow to give a shape to a cloud of points and are widely used in computational geometry. The use of αshapes in the context of NEM allows the construction of models entirely in terms of nodes and also ensures the linear precision of the shape function even over nonconvex boundaries. This is particularly true in the case of non-strongly non-convex domains. Indeed as mentioned by Cueto et al. [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF][START_REF] Cueto | Overview and recent advances in Natural Neighbour Galerkin Methods[END_REF], the α-shapes, used in such a method, fail to ensure the linearity of the shape functions over highly non-convex boundaries (Fig. 3.6), e. g. around a crack tip. α-NEM is therefore limited to nonstrongly concave geometries. Otherwise, the same problem related to the imposition of Dirichlet boundary conditions in the NEM, is raised. More recently, Yvonnet et al. [START_REF] Yvonnet | The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces[END_REF][START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous problems: the constrained natural element method (C-NEM)[END_REF][START_REF] Yvonnet | Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoi constraint-Méthode des éléments C-Naturels[END_REF], have proposed another extension of the NEM to ensure the linearity of the shape functions over all types of boundaries (even strongly non-convex): the Constrained Natural Element Method (CNEM). The CNEM can be regarded as another extension of the NEM which retains its attractive properties and overcomes the difficulties related to its application for strongly nonconvex domains. The main originality of this approach is the introduction of a visibility criterion [START_REF] Yvonnet | The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces[END_REF][START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous problems: the constrained natural element method (C-NEM)[END_REF][START_REF] Yvonnet | Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoi constraint-Méthode des éléments C-Naturels[END_REF][START_REF] Illoul | Mise en oeuvre de la méthode des éléments naturels contrainte en 3D Application au cisaillage adiabatique[END_REF][START_REF] Yvonnet | Nouvelles approches sans maillage basées sur la méthode des éléments naturels pour la simulation numérique des procédés de mise en forme[END_REF] in the NEM to select natural neighbors for constructing the interpolation. A modified so-called constrained Voronoï diagram is introduced for the computation of the shape functions. As a meshless method, CNEM has practically all the advantages of this class of methods. Also, with the introduction of a visibility criterion to select the natural neighbors in the construction of shape functions, it has overcome almost all the disadvantages of the meshless methods.

Classification and choice of continuum method: CNEM

The grid-based methods have been widely applied in various areas of CSM and CFD, and are the dominant methods in Computational Mechanics. Table 3.1 presents a comparison between the different categories of the grid-based methods.

As shown in this table, the Eulerian methods and the Combined Lagrangian-Eulerian methods represent the "weak links" in the grid-based methods with respect to the objectives of this dissertation. Indeed, the advantages of these methods are largely related to large and extremely large deformations. However, the discretecontinuum coupling method which will be developed in this work is mainly devoted to Computational Solid Mechanics where deformations are not so significant. Besides, the Discrete Element Method which is retained as a part of the coupling approach is based on the Lagrangian description. Hence, the coupling of this method with a continuum method using the same description seems to be advantageous.

The Lagrangian methods which are mainly represented by the Finite Element Method (FEM) are based on Lagrangian description and are more suited for CSM simulations (Tab. 3.1). Bearing in mind that the sought continuum method will be applied far from the regions of complex phenomena (which will be modeled by the DEM), the Finite Element Method represents a good candidate, among the grid-based methods, to be coupled with the DEM.

As for the meshless methods, the Lagrangian description is used in almost all the methods belonging to this class. Table 3

.2 gives a comparison between the different

Lagrangian methods

Eulerian methods Combined L-E Discretization

Mesh attached to material Mesh fixed in the space

Updated mesh

Examples FEM [START_REF] Zienkiewicz | Finite Element Method: Its Basis & Fundamentals[END_REF][START_REF] Zienkiewicz | The Finite Element Method for Solid and Structural Mechanics[END_REF][START_REF] Gouri | Une présentation de la méthode des éléments finis[END_REF] FVM [START_REF] Versteeg | An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Pearson education limited[END_REF] CEL [START_REF] Mair | Hydrocode methodologies for underwater explosion structure medium/interaction[END_REF][START_REF] Xie | Application of a coupled Eulerian-Lagrangian method to simulate interactions between deformable composite structures and compressible multiphase flow[END_REF][START_REF] Felici | A coupled Eulerian/Lagrangian method for the solution of threedimensinal vortical flows[END_REF]] X-FEM [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] FDM [ The main drawbacks of the Approximation methods are related to the treatment of the problem domain boundaries, especially the imposition of essential boundary conditions (e.g. displacements). The choice of an approximation method as a part of the coupling approach can, therefore, entail more difficulties in the coupling implementation. Indeed, to ensure correct communication between the coupled methods, it is important to correctly treat the interfaces between them. Hence, additional efforts must be undertaken when approximation method is used. Therefore, this class of methods will not be retained for the competition.

The Interpolation methods have overcome almost all the difficulties encountered in Approximation methods. The first method in this class which is the Natural Element Method (NEM) is mainly adapted for convex domains. Also, the ff-NEM, which is an extension of the NEM, is more suited for non-extremely non-convex domains. With the aim to develop a generic coupling method able to deal with all types of domains (even extremely non-convex), these two methods will not be retained. The CNEM which is also an extension of the NEM has circumvented the geometry limitations, while keeping the attractive properties of the factory method (NEM). Therefore, it is a good candidate to be coupled with the DEM.

Subclasses

Approximation methods

Interpolation methods

Discretization

Scattered nodes attached to the material

Scattered nodes attached to the material

Examples SPH [START_REF] Lucy | Numerical approach to testing the fission hypothesis[END_REF][START_REF] Liu | Smoothed particle hydrodynamics : a meshfree particle method[END_REF] NEM PUFEM [START_REF] Melenka | The partition of unity finite element method: Basic theory and applications[END_REF] α-NEM EFGM [START_REF] Belytschko | Element-Free Galerkin methods[END_REF] CNEM MLPG [START_REF] Atluri | A new meshless local Petrov-Galerkin (MPLG) approach in computational mechanics[END_REF] RKPM [START_REF] Li | Meshfree and particle methods and their applications[END_REF] PIM [START_REF] Liu | A point interpolation method for two-dimensional solids[END_REF] FPM [START_REF] Oñate | A Stabilized Finite Point Method for Analysis of Fluid Mechanics Problems[END_REF][START_REF] Oñate | A Mesh-Free Finite Point Method for Advective-Diffusive Transport and Fluid Flow Problems[END_REF][START_REF] Oñate | A Finite Point Method for Incompressible Flow Problems[END_REF] Then, at this stage, two numerical methods present good candidates to be coupled with the DEM: a grid-based method (FEM) and a meshless method (CNEM). Hereafter, a comparison between these methods will be given to select the most appropriate.

Track

Movement of material points

Movement of material points

Implementation

⋆⋆⋆ ⋆⋆

Table 3.3 compares the FEM and CNEM on some fundamental aspects. These two methods are very related through the following points : (i) the need of an underlying structure; (ii) the boundary conditions can be directly applied, due to the interpolation character and the linearity of the shape functions at the boundaries and interfaces. Applied on the same mechanical problem, the FEM computation is generally less expensive, in terms of CPU time, and more robust than the CNEM computation [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF]. Indeed, a step for construction of Delaunay tessellations and the computation of shape functions is needed in the case of CNEM which can be costly. Besides, the use of rational shape functions in CNEM (vs. polynomial shape functions in FEM) makes the CNEM integration more laborious and less robust. However, since only scattered nodes (without any predefined connectivity between them) are used in CNEM for spatial discretization, this method is more suited to discretize complex geometries. Also, CNEM is more advantageous to treat relatively large deformations characterizing the response of solids in highly dynamics, which is the scope of this work. Another point fostering the choice of the CNEM, compared to the FEM, is due to its "meshless" nature. Since, DEM is also a meshless method, the coupling between CNEM and DEM seems more natural. No additional efforts are required, due to the presence of mesh at the interface between these methods.

For these reasons, the CNEM is finally retained to be coupled with the DEM in this work. The main specificities of this method will be reviewed later in Chapter 4.

FEM CNEM Discretization

⋆⋆ ⋆⋆⋆

Implementation

⋆⋆⋆ ⋆⋆

Conclusion

As seen above, there are many numerical methods used in Computational Mechanics. These methods can mainly be classified into: discrete methods based on Newtonian mechanics; and continuum methods based on continuum mechanics. The first class is very time-consuming but more adapted to deal with complex problems such as fracture, wear and abrasion. The second class, which is predominant in both CSM and CFD, is less expensive in terms of CPU than the first one, however, it comes up against its limits when faced to problems with discontinuities. In this case, additional treatments are required to deal with these problems which can result in losing the advantages of this class of methods. The advantages of these two classes can largely be regarded as complementary. Therefore, it would be benifical to combine two methods, belonging to each of these classes: discrete method and continuum method.

The main objective of this chapter is to select the numerical methods that will be coupled in this work. To simplify the choice of the appropriate methods, the above classes were divided into several subclasses containing methods that share fundamental features. The discrete methods were divided into three subclasses: Lattice models, particle methods and contact dynamics (Fig. 3.1). Because the discrete-continuum coupling approach is attempted to models highly dynamic problems, the contact dynamics subclass, which is adapted for static and non-highly dynamic problems, was ignored. Then, an "hybrid" method combining the features of the Lattice and particle models was retained. Specifically, the variant of the Discrete Element Method recently developed by André et al. [4] is chosen to be coupled with a continuum method (unknown at this stage).

The continuum methods were divided into subclasses: grid-based methods and meshless methods (Fig. 3.2). The first subclass was divided into Lagrangian, Eulerian and Combined Lagrangian-Eulerian methods, according to the description frame of the governing equations. The two last categories are generally used for fluid dynamic problems which are not the scope of this dissertation. Also, since the chosen discrete method is based on the Lagrangian description, it is more natural to couple it to a Lagrangian continuum method. Therefore, only the Lagrangian grid-based methods were retained at this stage. The sought continuum method will be applied far from the high loaded regions, where complex phenomena can take place. Therefore, the FEM presents a good candidate for the grid-based methods.

The meshless methods were classified into Approximation methods and Interpolation methods. The first category suffers from several difficulties related to the treatment of the boundaries, which can make the coupling implementation more difficult. The Approximation methods were then ignored. From the Interpolation methods, the CNEM seemed to be the best approach to be coupled with the DEM. Indeed, it is the most recent interpolation method and has practically circumvented all the drawbacks of the meshless methods. In addition to the FEM, CNEM can also be a good candidate to be coupled with the DEM.

Although the FEM computation is generally more robust than CNEM computation, the advantages of the CNEM are mainly manifested in the fact that this method is more suited to discretize complex geometries and to treat relatively large deformations. Also, coupling DEM with CNEM would be more natural and easier than when using FEM. Because, both DEM and CNEM are meshless, no additional efforts are required due to the presence of mesh at the interface between the two methods. For these reasons, the CNEM was retained, as the continuum method, to be coupled with the DEM in this dissertation.

In conclusion, the two methods that were selected in this chapter are: the Discrete Element Method (DEM) from the discrete methods; and the Constrained Natural Element Method (CNEM) from the continuum methods. The main specificities of these methods as well as the developement of the coupling approach will be detailed in the next Chapter "Development of the discrete-continuum coupling method: DEM-CNEM".

Chapter 4

Development of the discrete-continuum coupling method: DEM-CNEM

"Science is a way of thinking much more than it is a body of knowledge."

Carl Sagan

Abstract

This chapter aims to develop a discrete-continuum coupling approach between the Discrete Element Method (DEM) and the Constrained Natural Element Method (CNEM), which have been selected in the previous chapter. First, the principle and the main specificities of these two methods will be recalled. Then, the coupling approach between these methods will be detailed. After, a parametric study of the different coupling parameters will be done. This study aims to simplify the setting of these parameters for practical applications. Finally, the new coupling approach will be validated using several reference dynamic tests. 57

Introduction

The simulation of the Laser Shock Processing on silica glass will be performed using a discrete-continuum coupling approach to be developed. Chapter 3 was devoted to the choice of the two methods that will be coupled. From the discrete methods, the variant of Discrete Element Method (DEM) recently implemented by André et al. [4,3] has been chosen. As for the continuum method, the Constrained Natural Element Method (CNEM) developed by Yvonnet et al. [START_REF] Yvonnet | A new extension of the natural element method for non-convex and discontinuous problems: the constrained natural element method (C-NEM)[END_REF][START_REF] Yvonnet | Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoi constraint-Méthode des éléments C-Naturels[END_REF][START_REF] Yvonnet | Nouvelles approches sans maillage basées sur la méthode des éléments naturels pour la simulation numérique des procédés de mise en forme[END_REF] has been retained. These two methods rely on different physical bases and have different scales. This chapter, which aims to couple these approaches, should therefore face two major challenges:

• The first challenge is due to the different physical bases of these methods. The DEM and CNEM variables do not have the same nature. Indeed, the discrete variables associated with the DEM are only defined in the particle positions. However, continuous field variables are sought in the case of CNEM. Therefore, a special treatment is required at the interface between the coupled methods to take on this challenge.

• The second challenge which mainly encountered in dynamics is due to the fact that these methods have different scales. In a multi-scale coupling approach, the high frequency portion of the propagating waves is often spuriously reflected at the small/coarse scale interface. This phenomenon has already been addressed using the Finite Element Method (FEM) with different element sizes [START_REF] Celep | Spurious reflection of elastic waves due to gradually changing finite element size[END_REF].

The importance of the multi-scale approaches has attracted many researchers. Therefore, numerous works have been published on the subject, and several coupling methods have been developed. Hereafter, some interesting coupling approaches in the literature will be briefly reviewed.

Ben Dhia [START_REF] Dhia | Analyse mathématique de la méthode Arlequin mixte[END_REF][START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Dhia | Problèmes mécanique multi-échelles: la méthode Arlequin[END_REF], in a pioneer work, has developed the Arlequin approach as a general framework which allows the intermixing of various mechanical models for structural analysis and computation. Abraham et al. [START_REF] Abraham | Spanning the continuum to quantum length scales in a dynamic simulation of brittle fracture[END_REF][START_REF] Broughton | Concurrent coupling of length scales: methodology and application[END_REF] developed a methodology that couples the tight-bending quantum mechanics with Molecular Dynamics (MD) such that the two Hamiltonians are averaged in an overlapping zone. A damping was used in the overlapping zone to reduce the spurious reflections at the interface between the two models. Nevertheless, the choice of the damping coefficient remains difficult. Smirnova et al. [START_REF] Smirnova | A combined molecular dynamics and finite element method technique applied to laser induced pressure wave propagation[END_REF] developed a combined Molecular Dynamics (MD) and Finite Element Method (FEM) model with a transition zone in which the FEM nodes coincide with the positions of the particles in the MD region. The particles in the transition zone interact with the MD region via the interaction potential. At the same time, they experience the nodal forces due to the FEM grid. Belytschko and Xiao [START_REF] Belytschko | Coupling methods for continuum model with molecular model[END_REF][START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF] have developed a coupling method for Molecular Dynamics and Continuum Mechanics models based on a bridging domain method. In this method, the two models are overlaid at the interface and constrained with a Lagrange multiplier model in the overlapping subdomain. Fish et al. [START_REF] Fish | Concurrent AtC coupling based on a blend of the continuum stress and the atomistic force[END_REF] formulated an atomistic-continuum coupling method based on a blend of the continuum stress and the atomistic force. In terms of equations, this method is very similar to the Arlequin method. In an interesting work, Chamoin et al. [START_REF] Chamoin | Ghost forces and spurious effects in atomic-to-continuum coupling methods by the Arlequin approach[END_REF] have analyzed the main spurious effects in the atomic-to-continuum coupling approaches and they proposed a corrective method based on the computation and injection of dead forces in the Arlequin formulation to offset these effects. Aubertin et al. [START_REF] Aubertin | A coupled molecular dynamics and extended finite element method for dynamic crack propagation[END_REF] applied the Arlequin approach to couple the extended finite element method X-FEM with the Molecular Dynamics MD to study dynamic crack propagation. Bauman et al. [START_REF] Bauman | Adaptive multiscale modeling of polymeric materials with Arlequin coupling and Goals algorithms[END_REF] developed a 3D multi-scale method, based on the Arlequin approach, between highly heterogeneous particle models and nonlinear elastic continuum models. Recently, Combescure et al. [START_REF] Chuzel-Marmot | Three dimentional SPH-FEM gluing for simulation of fast impacts on concrete slabs[END_REF] formulated a 3D coupling method, applied for fast transient simulations, between the SPH and the FEM. This coupling method is, also, based on the Arlequin approach. For more details, a review of these methods can be found in Reference [START_REF] Curtin | Atomistic /continuum coupling in computational materials science[END_REF].

The different coupling approaches can be divided in two classes: edge-to-edge methods and methods with overlapping zones, called overlapping methods (Fig.

2.25).

The first class [START_REF] Belytschko | Coupling methods for continuum model with molecular model[END_REF][START_REF] Tadmor | Quasicontinuum analysis of defects in solids[END_REF] is mainly applied for static studies. Indeed, using an edge-to-edge method, it is very difficult to reduce spurious reflections at the interface between models. The second class seems to be more applicable to dynamic studies, which is the scope of this dissertation. Therefore, a discrete-continuum coupling approach with overlapping region is developed in this work. This coupling approach is based on the Arlequin formulation which allows high frequencies damping and avoids spurious wave reflections without any additional filtering or artificial damping. Indeed, the fine scale solution (of the discrete model) is projected onto the coarse scale solution (of the continuum model) in the overlapping zone at each time step. This filters the high frequency waves coming from the fine scale model (discrete model), which are larger than the cutoff frequency of the coarse scale model (continuum model).

After recalling the main specificities and features of the DEM and CNEM, this chapter will detail how these methods are coupled. Then, the different coupling parameters involved in this coupling approach will be studied to simplify its application in practice. Finally, the DEM-CNEM coupling method will be validated using several reference dynamic tests: dynamic tension, bending, torsion and impact tests on 3D beams.

The discrete element method: DEM

The variant of the DEM used in this work is that implemented by André et al. [3,4]. It can be regarded as an "hybrid" method between Lattice and Particle models. In this method, a continuum is represented by a set of particles (discrete elements) having volumes and masses. 3D cohesive beam bonds are used to link these particles. These beam bonds allow the transmission of tensile, bending and torsional efforts to the particles. Figure 4.1 presents an illustration of DEM modeling of a continuum. André et al. [3,4] have shown that it is possible to obtain under certain conditions the equivalence of the geometric and mechanical properties between a continuum and its DEM representation. These conditions can be divided in two classes: conditions on the geometric discretization and conditions on the mechanical properties of

Construction of the DEM domain

The construction of a DEM domain must take into account the structural properties of the problem domain. For the case of silica glass, two main structural properties must be considered: homogeneity and isotropy. Therefore, the DEM representation of this material must be validated with regard to these two properties.

Construction technique

It is shown in the literature that the DEM representation of a domain can influence its mechanical state. For example, the use of an orderly geometric configuration can lead to anisotropic mechanical behavior [START_REF] Schlangen | New method for simulating fracture using an elastically uniform random geometry lattice[END_REF][START_REF] Carmona | Fragmentation processes in impact of spheres[END_REF]. Since silica glass is an isotropic material, disorderly configurations must be used for the DEM representation of this material. Different techniques can be found in the literature, used to construct DEM domains. Among them, one can cite the dynamic compaction [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF][START_REF] Martin | Study of particle rearrangement during powder compaction by the discrete element method[END_REF]. Although this technique is widely used in the literature, it is restricted to relatively simple domains. In this work, we have used another technique which consists of three steps:

• Random free filling: the boundaries of the DEM domain must first be drawn, which are the same as in the real domain. Subsequently, discrete elements are randomly placed inside the hull formed by these boundaries, such that no interpenetration is allowed between them. To obtain disorderly DEM domains, a statistical distribution X (generally uniform distribution) is applied on the discrete elements radii. This distribution is defined as:

X = r max -r min r mean (4.1)
where r max , r min and r mean are the maximum, minimum and mean radii of the particles. This first step achieves when all the free places become occupied by discrete elements (i.e. it becomes impossible to place a new discrete element without interpenetration). Figure 4.2 illustrates the evolution of the DEM domain filling during this first step.

• Tamping: this step aims to continue the DEM domain filling so that to obtain a fully compacted isotropic domain. An iterative process is run to randomly place new discrete elements. At each iteration, a discrete element is randomly placed within the pre-filled DEM domain; then, a DEM computation is performed to re-organize the discrete elements and minimize the interpenetrations between them. This step is very time-consuming since the DEM computation must continue until a stabilization criterion is fulfilled (the total kinetic energy must be smaller than a certain limit). A solution to reduce the computation time is to place a set of n discrete elements at a time. Research is currently underway to determine the optimal value of n.

• Relaxation: at the end of the second step, the DEM domain has a relatively high mean interpenetration δ between the discrete elements. In the case of granular media simulations where the discrete element forces and torques are determined only by contacts, it is important to reduce this parameter (δ). This is performed by reducing gradually the stiffness of the confining hull. In this work, since the computation of the discrete element forces and torques are based on the cohesive beams (and not on the contact laws between discrete elements), this step is not necessary.

Homogeneity According to previous works in the literature [START_REF] Gotoh | Statistical geometrical approach to random packing density of equal spheres[END_REF][START_REF] Finney | Random packings and the structure of simple liquids. i. the geometry of randomclose packing[END_REF], to satisfy the homogeneity property of the constructed domain, two criteria must be fulfilled. The first one concerns the coordination number n coord (average number of discrete element neighbors) which must be around 6 [START_REF] Gotoh | Statistical geometrical approach to random packing density of equal spheres[END_REF]; the second one concerns the volume fraction, noted v f , which must be close to 0.63 [START_REF] Finney | Random packings and the structure of simple liquids. i. the geometry of randomclose packing[END_REF]. André et al. [3,4] have shown that using a uniform distribution X = 25% on the particles radii, highly irregular DEM domain can be obtain and these criteria can be fulfilled. Figure 4.3 presents two DEM geometric configurations obtained for X = 0% and X = 25%.

Isotropy Another geometric property which must be verified to ensure correct mechanical behavior is the geometric anisotropy. This property is defined as the repartition of the contact directions in the space. If this repartition is perfectly uniform, the associated material can be considered as isotropic. To measure the geometric anisotropy, André et al. [3,4] have proposed a geometric approach using a platonic solid (geode) which discretizes the 4 π sr (steradian) full 3D space into 320 equal facets. Each facet corresponds to an elementary solid angle of 4 π/320 sr (Fig. 4.4). This geode is used as the support of a 3D histogram for the contact directions in the DEM domain.

Each contact between two discrete elements i and j is characterized by the directionvector O i O j that joins the discrete element centers. The contact line passing through the geode center and parallel to O i O j intersects two opposite facets. The considered contact is then attached to these facets. Repeating this process for all the contacts in the DEM domain allows the computation of the facets weights (number of contacts attached to these facets). The geometric anisotropy is characterized by the standard deviation of these weights. Figure 4.5 presents the 3D histograms of the contact directions weights for X = 0% and X = 25%. Note that in the case of a perfectly uniform distribution of the contact directions, the associated 3D histogram has the shape of a geode (Fig. 4.4). For X = 0%, the compacted domain is highly anisotropic. However, the distribution X = 25% provides a better representation for an isotropic domain, and then it will be used in this work to construct the silica glass numerical samples.

DEM domain fineness

The fineness of a DEM domain is characterized by the number of particles (discrete elements) n p which are used to represent the geometrical domain. Three criteria are used to drive the homogeneity and isotropy of the discrete domain (DEM domain): the coordination number n coord , the volume fraction v f and the standard deviation of the weights of the contact directions. The question that arises here is what is the minimum n p which must be used to ensure stable geometric properties. According to the works of André et al. [4], beyond n p = 10 000,

Cohesive beam bond model

After construction of the DEM domain, the cohesive links are introduced between the neighbor particles to ensure the cohesion of the medium and to model its behav- ior. To distinguish the microscopic properties of these links from the macroscopic properties of the full problem, the microscopic parameters will be subscripted by µ.

In this DEM variant, cylindrical cohesive beams are used to link the neighbor discrete elements. Therefore, only two independent geometric parameters are required to completely define the beams geometry: r µ and l µ (radius and length of the beams, respectively). To characterize the mechanical behavior of these beams, two mechanical parameters are used: a microscopic Young's modulus E µ and a microscopic Poisson's ratio ν µ . No mass properties are assigned to the beams, these properties are only assigned to the discrete elements. Figure 4.6 presents an illustration of a cohesive beam between two discrete elements. The Euler-Bernoulli analytical beam theory is used to compute the beam forces and torques on the discrete elements connected to this beam ((4.2), (4.3), (4.4) and (4.5)).

F DE1 µ = +E µ S µ ∆l µ l µ x + 6E µ I µ l 2 µ (-(θ 2z + θ 1z )y + (θ 2y + θ 1y )z) (4.2) F DE2 µ = -E µ S µ ∆l µ l µ x - 6E µ I µ l 2 µ (-(θ 2z + θ 1z )y + (θ 2y + θ 1y )z) (4.3) T DE1 µ = + G µ Io µ l µ (θ 2x + θ 1x )x - 2E µ I µ l µ ((θ 2y + 2θ 1y )y + (θ 2z + 2θ 1z )z) (4.4) T DE2 µ = - G µ Io µ l µ (θ 2x + θ 1x )x - 2E µ I µ l µ ((2θ 2y + θ 1y )y + (2θ 2z + θ 1z )z) (4.5)
Where:

• R(x, y, z) is the beam local frame where x is the beam axis. Table 4.1: The silica glass microscopic properties of the cohesive beam bonds; rµ is a dimensionless cohesive beam radius, defined as the ratio between the beam radius and the mean particle radius.

• F DE1 µ and F DE2 µ are the beam force reactions acting on discrete elements 1 and 2 (connected to this beam bond).

• T DE1 µ and T DE2 µ are the beam torque reactions acting on discrete elements 1 and 2.

• l µ and ∆l µ are the initial beam length and the longitudinal extension.

• θ 1 (θ 1x , θ 1y , θ 1z ) and θ 2 (θ 2x , θ 2y , θ 2z ) are the rotations of the beam cross sections expressed in the beam local frame.

• S µ , Io µ and I µ are respectively the beam cross-sectional area, the polar moment of inertia and the second moment of area with respect to the y and z axes.

• E µ and G µ are respectively the Young's and shear moduli.

Table 4.1 presents the microscopic properties of the elastic cohesive beams which ensure the expected macroscopic elastic behavior of silica glass. These properties have been determined by numerical calibration tests as explained in Reference [4].

The constrained natural element method: CNEM

As an extension of the NEM, the Constrained Natural Element Method (CNEM) is very similar to the factory method (NEM) and keeps all its attractive properties. This method differs only in the introduction of a visibility criterion to restrict (select) the natural neighbors. This has allowed to circumvent the problem of the imposition of essential boundary conditions over all types of boundaries and interfaces (convex or not).

For clarity, the main features of the factory method (NEM) will first be reviewed. Then, the visibility criterion as well as the resulting modifications will be introduced.

Natural Neighbor (NN) interpolation

The notion of natural neighbors and natural neighbor (NN) interpolation were introduced by Sibson [161] for data fitting and smoothing of highly irregular set of nodes. The NN interpolation is based on the Voronoï diagram and the Delaunay tessellations which were originally proposed by mathematicians [START_REF] Delaunay | Sur la sphère vide. A la mémoire de Georges Voronoï[END_REF][START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques[END_REF] and later applied in several engineering fields.

Voronoï diagram

The Voronoï diagram is defined as a set of cells T i associated to the nodes i such that all points in a particular cell are closer to its defining node than any other one (Fig.

4.7).

The Voronoï diagram of a given set of nodes is unique, it fills the convex hull of these nodes. Mathematically, this diagram is defined in any dimension d as follows:

T i = x ∈ R d , xx i < xx j , ∀ i = j , i = 1..n t (4.6)
where n t is the total number of nodes and xx i represents the distance between a point of x coordinates and the node i located in x i .

In duality with the Voronoï diagram, the concept of the Delaunay tessellation (or triangulation in 2D) can be introduced. This concept was proposed by Voronoï [START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques[END_REF] and extended by Delaunay [START_REF] Delaunay | Sur la sphère vide. A la mémoire de Georges Voronoï[END_REF]. It consists in connecting the nodes whose the Voronoï cells share some boundaries to construct the Delaunay tessellations (triangles in 2D, tetrahedrons in 3D). An important property of the Delaunay tessellations is that the Circumscribed spheres (mathematical definition of a sphere is used here, this is, circle in 2D or sphere in 3D) of these elements (Delaunay tessellations) do not contain any node. The natural neighbors of a node i are those connected to i by an edge in the Delaunay tessellations. In other words, nodes having Voronoï cells which share some boundaries are natural neighbors. If the nodes are relatively spaced in some regions or distributed in highly anisotropic manner, the set of natural neighbors will also reflect these features. However, it always represents the best set of surrounding neighbors. Natural neighbors are natural candidates to define a basis of shape function construction and local interpolation scheme.

NN shape functions

To define the neighbor relation of any point placed into the Delaunay tessellation, Sibson [START_REF] Sibson | A vector Identity for the Dirichlet tessellations[END_REF] used the concept of second-order Voronoï diagram which is defined as a set of cells T ij . A cell T ij is the locus of the points x having i as the closest node and j as the second closest node. Mathematically, this diagram is defined in any

T ij = x ∈ R d , xx i < xx j < xx k , ∀ i = j = k , i = 1..n t (4.7)
where n t is the total number of nodes and xx i represents the distance between a point of x coordinates and the node i located in x i . Several types of shape functions are used in the NN interpolation.

Sibson shape functions

Assuming that K denotes the Lesbegue measure (length, area or volume in respectively 1D, 2D or 3D), K x denotes the Lesbegue measure of T x and K xi denotes the Lesbegue measure of T xi . The Sibson shape function φ i associated to a node i is defined as the ratio of K xi and K x , where x are the coordinates of the point at which the shape function is calculated.

φ i (x) = K xi K x (4.8)
In (4.8), K x can also be calculated as

K x = n v ∑ i=1 K xi ,
φ i (x) = area A area B (4.9)
The Sibson shape functions have C 0 continuity at the nodes, C 1 continuity at the Delaunay spheres (circles in 2D) and C ∞ continuity anywhere else.

Laplace (non-sibsonian) shape functions Others NN shape functions, called Laplace or non-sibsonian shape functions, are used in the NN interpolation [START_REF] Belikov | The non-Sibsonian interpolation : a new method of interpolation of the values of a function on an arbitrary set of points[END_REF]. Assuming that f xi is the Lesbegue measure of the Voronoï facet (length in 2D or area in 3D) separating the point of x coordinates and the node i; and d xi is the distance between them. The Laplace shape function associated to i is calculated in x as:

φ i (x) = α xi n v ∑ i=1 α xj , with α xi = f xi d xi (4.10)
where n v is the number of the natural neighbors of x. The right view of Figure 4.8 gives and illustration of 2D calculation of Laplace shape functions. The Laplace shape functions have C 0 continuity at the nodes and at the Delaunay spheres (circles in 2D), and C ∞ continuity anywhere else. These functions are less expensive in terms of CPU than the Sibson shape functions. Because only areas (lengths in 2D) require to be evaluated in these shape functions, whereas, volumes (area in 2D) must be calculated in the Sibson shape functions. However, the Sibson shape functions are more regular at the Delaunay spheres (circles in 2D) (C 1 continuity).

Other NN shape functions A generalization of NN-based shape functions, called k thorder standard coordinates, has been proposed by Hiyoshi and Sugihara [START_REF] Hiyoshi | Improving continuity of Voronoi-based interpolation over Delaunay spheres[END_REF]. These functions have C 0 continuity at the nodes, C k continuity (k is a nonnegative integer) at the Delaunay spheres and C ∞ continuity anywhere else. The cases k = 0 and k = 1 coincide respectively with the Laplace and Sibson shape functions. It should be noted that the computation time of the shape functions is proportional to the continuity order k at the Delaunay spheres. Therefore, using highly regularized shape functions (large k) is very time-consuming and can become crippling.

Remark In this work, Sibson shape functions have been retained for interpolation. These functions are more regular than the Laplace shape functions. Also, they are not very costly in terms of CPU compared to the k thorder standard coordinates with (k > 1). Based on the above shape functions, the NN interpolation of a vectorvalued function u can be derived:

u h (x) = n v ∑ i=1 φ i (x) u i (4.11)
where n v is the number of natural neighbor nodes of the point having x coordinates, φ i is the Sibson shape function associated to a node i and u i is the nodal value of u at the node i of x i coordinates (u i = u(x i )).

Support of NN shape functions

The support of a NN shape function φ i (x) associated to i is defined as the space occupied by the union of the n v circumspheres (circumcircles in 2D) of the Delaunay tetrahedrons (triangles in 2D) in which i is a vertex; n v is the number of the natural neighbors of i. represented by the best set of surrounding neighbors even in the case of irregular of highly anisotropic distribution of nodes.

Properties of NN shape functions

The main properties of the NN shape functions are:

• The Kronecker delta property: This property states that the NN interpolation passes through the nodal values.

φ i (x j ) = 1 i = j 0 i = j (4.12)
• Compact support: the NN shape functions are compactly supported. It means that they are nonzero only in delimited regions and vanish beyond. The NN interpolation is therefore transformed from global operation to local operation. This will lead to a set of sparse discretized system matrices, which considerably reduces the computational efforts.

• Partition of unity: by construction, the NN shape functions verify the following relation:

n t ∑ i=1 φ i (x) = 1 ∀ x ∈ Ω (4.13)
where n t is the total number of nodes. This property is also called 0 thorder (C 0 ) consistency and means that the NN interpolation can exactly reproduce constant functions.

• Linear consistency: the NN shape functions also satisfy the local coordinate property (1 thorder (C 1 ) consistency) defined as:

x = n t ∑ i=1 φ i (x)x i ∀ x ∈ Ω (4.14)
This property combined with the partition of unity property implies that the NN interpolation have the linear completeness property. The NN interpolation exactly reproduces any linear polynomial. It should be noted that in solving of PDEs based on a weak formulation (e.g. the Galerkin formulation), there is a minimum consistency requirement to ensure the convergence of the numerical results, depending on PDEs order. For the Galerkin formulation, C k consistency is required to solve PDEs of 2k order.

• Linear variation on convex boundaries: the NN shape functions are stricly linear on convex boundaries. Combined with the Kronecker delta property, this property implies that the essential boundary conditions can directly be applied. However, in the case of non-convex boundaries, this property does not satisfied. This problem is circumvented by introducing the notion of visibility criterion and then the Constrained Natural Neighbor (CNN) interpolation.

Visibility criterion

This criterion states that a point of x coordinates is in the influence domain of a node i if x is within the region where the shape function φ i associated to i is non-null (φ i (x) = 0) and it is visible from the node i when the boundaries are assumed opaque. Assuming that Ω C is the whole continuum domain and

Ω i = {x|x ∈ Ω C , φ i (x) = 0}
is the original (convex) influence domain of a node i, this criterion is defined mathematically as:

Ω visible i = {x|x ∈ Ω i , l x→x i ∩ ∂Ω i = ∅} (4.15)
where l x→x i represents the straight line relying x and x i (coordinates of i) and ∂Ω i is the boundary of Ω i .

Constrained Natural Neighbor (CNN) interpolation

Taking in account the visibility criterion, the Voronoï diagram is extended to become the constrained Voronoï diagram [START_REF] Edelsbrunner | Geometry and Topology for Mesh Generation[END_REF] which is defined as:

T C i = x ∈ R d , xx i < xx j , ∀ i = j, l x→x i ∩ ∂Ω i = ∅, l x→x j ∩ ∂Ω j = ∅ , i = 1..n t (4.16)
Concerning the Constrained Delaunay tessellation, it is always possible to construct the constrained Delaunay triangulation for a given set of nodes in the plane (in 2D). However, as reported by several researchers [START_REF] Shewchuck | Tetrahedral mesh generation by Delaunay refinement[END_REF][START_REF] Yvonnet | Nouvelles approches sans maillage basées sur la méthode des éléments naturels pour la simulation numérique des procédés de mise en forme[END_REF], its generalization to 3D does not always exist without adding new nodes. Some techniques have been proposed in the literature to construct the constrained Delaunay tessellations in 3D by addition of Stainer points [START_REF] Shewchuck | Tetrahedral mesh generation by Delaunay refinement[END_REF].

With introduction of the visibility criterion, the natural neighbors of any point are confined to the constrained natural neighbors which are visible by the considered point. Once the constrained Voronoï diagram as well as the associated Delaunay tessellation are constructed, the selection of the constrained natural neighbors can directly be determined, following the same procedure that in NEM. The associated CNN shape functions φ C i can also be determined using the constrained Voronoï diagram and the constrained second-order Voronoï diagram. In such diagrams, any Voronoï cell facet related to two non-mutually-visible nodes is removed. In this case, the CNN shape functions φ C i restore the linearity property on any type of boundaries. The NN interpolation given by (4.11) becomes:

u h (x) = n C v ∑ i=1 φ C i (x) u i (4.17)
where n C v is the number of the CNN nodes.

Numerical integration

In Computational Mechanics, the global weak formulation of the studied problem is derived, based on several approaches such as the Galerkin approach (used in CNEM). This formulation is given in continuous form which is very difficult, if not possible, to be analytically solved. However, approximated solutions can be obtained numerically by changing the continuous weak formulation into a discrete system of equations. The integrals in this system are evaluated numerically using several techniques. Among them, one can cite the Gauss Integration and the Direct Nodal Integration. In the case of CNEM, the Gauss technique might be inaccurate since the shape functions are not generally polynomials and the supports of these functions do not coincide with the integration cells [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF][START_REF] González | Numerical integration in Natural Neighbour Galarkin methods[END_REF]. Also, direct nodal integration , using the nodes as integration points, typically leads to numerical instabilities resulting from the under-integration of the weak formulation as reported in the literature [START_REF] Beissel | Nodal integration of the element-free Galerkin method[END_REF]. Recently, new numerical integration procedures have been proposed in the framework of meshless methods. In this work, the Stabilized Conforming Nodal Integration (SCNI) developed by Chen et al. [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF] will be used to compute the integral terms coming from the continuum model. This technique which is origi-nally developed for strain localization problems relies on the gradient smoothing to stabilize the nodal integration as follows:

∇u h (x i ) = ˆΩi w(x, x -x i ) ∇u h (x i ) dΩ (4.18)
where x i are the coordinates of the node i, Ω i is a subdomain surrounding i (generally taken as the Voronoï cell associated to i). w is a weight function having a compact support (vanishes beyond a certain region) and satisfying these properties: w ≥ 0 and ´Ωi w(x, xx i ) dΩ = 1. This function is generally taken as the step function defined on Ω i , supposing that V i is the Lesbegue measure of Ω i :

w(x, x -x i ) = 1 V i if x ∈ Ω i 0 if x / ∈ Ω i (4.19)
Using (4.18) and (4.19), the smoothed strain at the node i can be determined:

εh kl (x i ) = 1 V i ˆΩi ε h kl (x) dΩ (4.20)
Applying the divergence theorem, the integral (4.20) can be transformed to hyperplane contour integral (the hyperplane is defined here as a subset of dimension d -1, d is the dimension of Ω i ):

εh kl (x i ) = 1 2V i ˆ∂Ω i (u h k n l + u h l n k ) dΓ (4.21)
n is the outward-pointing unit normal vector to ∂Ω i . The contour integrals are then approximated by successive application of Gauss Integration technique to each facet of Ω i . Introducing now the CNN interpolation, i.e. replacing u h by its expression given by (4.17), the strain tensor at i, expressed in engineering notation, can be written: εh is the smoothed gradient matrix [START_REF] Chen | A stabilized conforming nodal integration for Galerkin mesh-free methods[END_REF][START_REF] González | Numerical integration in Natural Neighbour Galarkin methods[END_REF] associated to l:

(x i ) = ∑ l∈Neighbors(i) [B l (x i )] {u l } (4.
[B l (x i )] =         b l1 (x i ) 0 0 0 b l2 (x i ) 0 0 0 b l3 (x i ) 0 b l3 (x i ) b l2 (x i ) b l3 (x i ) 0 b l1 (x i ) b l2 (x i ) b l1 (x i ) 0         , b lk (x i ) = 1 V i ˆ∂Ω i φ C l (x) n k dΓ (4.23)

Discrete-continuum coupling method: DEM-CNEM

In the previous sections ( §4.2 and §4.3), the main features of the retained methods have been reviewed. The present section will detail how these methods are coupled.

Arlequin approach: brief description

As discussed earlier, the coupling approach developed in this work is of overlapping type (i.e. the two methods are overlaid at the interface). This approach is based on the Arlequin technique which gives a flexible framework to couple dissimilar methods. This technique is mainly made up of three steps:

1. Decomposition of the global problem domain into two (or more) subdomains with overlapping regions; and superposition of the mechanical states in the resulting subdomains (Fig. 4.10).

A weak coupling (based on the weak formulation): the second point of the

Arlequin approach states that, to couple the superposed models, the discrepancy between the mechanical states, e.g. displacement, deformation, strain, etc, must be controlled using some kinds of fictive forces. However, to allow each model to express its own wealth, this discrepancy should be controlled in a weak manner and using averaging operators. This point consists in four steps: a) Definition of the gluing zone (Ω OG ): in a generic Arlequin coupling, the control of the mechanical states can be made only in a portion of the overlapping region (Fig. 4.10). This portion is named gluing zone.

b) Mediator space: only certain elements serve to compare the mechanical states in the gluing zone. These elements form a physical space, called "mediator space" (M). d) The junction model: the control of the mechanical states discrepancy in the gluing zone is made using a so-called junction model which exerts coupling forces between the coexistent models in this region. The junction models can be of different natures : rigid models such as the Lagrangian model, elastic models such as the penalty model or also a combination of these models such as the penalty-duality model.

3. Energy partition between the superposed models: as shown in Figure 4.10, the coupled models coexist in the overlapping region. Therefore, the energies in this zone must be weighted (multiplied) by a kind of partition of unity functions, called weight functions.

Arlequin approach: application to the DEM-CNEM coupling

A domain Ω G is considered with boundary ∂Ω = ∂Ω u + ∂Ω t such that essential (u g ) and natural (t g ) boundary conditions are prescribed respectively on ∂Ω u and ∂Ω t . For the sake of clarity, Ω G is supposed to be a 2D rectangular domain, however, all the following analyses are valid for any other geometry. This domain is divided into two subdomains, Ω C and Ω D , with an overlapping zone. These subdomains are modeled using respectively the continuum (CNEM) and discrete (DEM) methods (Fig. 4.12). An isotropic linear elastic behavior and small strain gradients are assumed for simplicity. The governing equations of both Ω D and Ω C are derived in Sections 4.4.2.1 and 4.4.2.2, while ignoring the coupling conditions. These conditions will be introduced in Section 4.4.2.3.

DEM formulation

In an isolated system of the discrete subdomain Ω D which is a set of n p spherical particles interacting via elastic cohesive beams, the governing equations can be written as: 3 such that, given the initial conditions:

i = 1..n p and t ∈ [0, t end ], find d i , θ i , f int i , c int i ∈ R 3 × R 3 × R 3 × R
f ext i + f int i = m i di c ext i + c int i = I i θi (4.24)
where d i , θ i , m i and I i represent respectively the displacement vector, the rotation vector, the mass and the mass moment of inertia of the particle i. f int i and c int i are respectively the internal force and torque applied on the particle i by its neighbors via the cohesive links. f ext i and c ext i are respectively the external force and torque applied on the particle i.

In keeping with the weak formulation in the continuum Ω C , Equation (4.24) can be reformulated as follows: Find (d, θ, f , c) ∈ D ad × O ad × F ad × C ad such that ∀ δ ḋ, δ θ ∈ Ḋad,0 × Ȯad,0 :

n p ∑ i=1 f ext i .δ ḋi + n p ∑ i=1 f int i .δ ḋi + n p ∑ i=1 c ext i .δ θi + n p ∑ i=1 c int i .δ θi = n p ∑ i=1 m i di .δ ḋi + n p ∑ i=1 I i θi .δ θi (4.25)
Where the admissible spaces are defined by:

• D ad = d = d i (t), i = 1..n p ∀ t ∈ [0, t end ] • O ad = θ = θ i (t), i = 1..n p ∀ t ∈ [0, t end ] • F ad = f = f int i (t), i = 1..n p ∀ t ∈ [0, t end ] • C ad = c = c int i (t), i = 1..n p ∀ t ∈ [0, t end ] • Ḋad,0 = ḋ = ḋi (t), i = 1..n p ; ḋ = 0 on ∂Ω u D ; ∀ t ∈ [0, t end ] • Ȯad,0 = θ = θi (t), i = 1..n p ; θ = 0 on ∂Ω u D ; ∀ t ∈ [0, t end ]

CNEM formulation

As an isolated system, the governing equations in the continuum subdomain Ω C can be written as:

∀ x ∈ Ω C (t) and t ∈ [0, t end ], find (u, σ) ∈ H 1 (Ω C ) 3 × L 2 (Ω C ) 6
such that, given the initial conditions:

           div(σ) + ρ f = ρ ü in Ω C σ = A : ε(u) ε(u) = 1 2 (∇u + ∇ t u) u = u g on ∂Ω u C σ.n = t g on ∂Ω t C (4.26)
where ρ is the density, u is the continuum displacement vector, σ is the Cauchy stress tensor, ε is the stain tensor, A is the stiffness tensor, f is the body force vector, u g and t g are respectively the given (prescribed) displacement and traction vectors on ∂Ω u C and ∂Ω t C , and n is the outward-pointing unit normal vector. The continuum weak formulation associated with Equation (4.26) can be written as follows: Find u ∈ U ad such that ∀ δ u ∈ U ad,0 , given the initial conditions:

ˆ∂Ω T C δ u • t g dΓ - ˆΩC ε(δ u) : A : ε(u) dΩ + ˆΩC ρ δ u • f dΩ = ˆΩC ρ δ u • ü dΩ (4.27)
where δ u represent the test functions and the admissible solution spaces, U ad and U ad,0 , are defined as follows:

• U ad = u = u(x, t) ∈ [H 1 (Ω C )] 3 ; u = u g on ∂Ω u C ; ∀ t ∈ [0, t end ] • U ad,0 = u = u(x, t) ∈ [H 1 (Ω C )] 3 ; u = 0 on ∂Ω u C ; ∀ t ∈ [0, t end ]

Coupling formulation

First, the gluing zone Ω OG where the mechanical states will be controlled must be defined. In this work, this zone is the same that the overlapping region Ω O : Ω OG = Ω O . Therefore, the control of the mechanical states will be made in whole Ω O . Hereafter, the term "overlapping region" will be used to design both Ω OG and Ω O .

To ensure a correct dialog between the discrete and continuum models, the control quantities in the overlapping zone must be carefully chosen. Regarding the present work, the velocities are chosen to be controlled, in a weak sense, in Ω O . As will be seen in Section 4.4.5, the velocity coupling (in which velocities are controlled in Ω O ) is easier and less expensive in terms of CPU than the displacement coupling (in which displacements are controlled in Ω O ). The mediator space M is therefore defined as the space of the velocities in Ω O .

To project the discrete and continuum velocities on the mediator space M, the CNN interpolation is used as projection operator Π. It should be recalled that Ω D is naturally a discrete system made up of a set of particles which are delimited by discontinuities. In other words, only discrete velocities at the particle positions can be defined. Therefore, interpolation of these quantities, in the strict sense, makes no physical meaning. To overcome this limitation, the part of the discrete domain restricted to Ω O (Ω D|Ω O ) will be considered as a continuum.

In the literature, several junction models are proposed in the Arlequin approach to couple dissimilar methods. In this work, the H 1 (Ω O ) scalar product is used to apply the coupling forces in the overlapping region. This model which is one of the most general junction models is defined as follows:

< λ, u -ḋ > H 1 (Ω O ) = ˆΩO λ • (Π u -Π ḋ) + l 2 ε(λ) : ε(Π u -Π ḋ) dΩ (4.28)
where (Π u -Π ḋ) is the difference between the projected continuum and discrete velocities on the mediator space M, λ is the Lagrange multiplier field and l is a parameter of the H 1 (Ω O ) junction model, named "junction parameter" in this manuscript. This parameter which has the dimension of a length is added to ensure the homogeneity of the integral terms in (4.28). In this work, l is considered as a variable which will be studied in Section 4.5. If l = 0, the H 1 (Ω O ) scalar product is reduced to the L 2 (Ω O ) scalar product which is well-known as the Lagrange multiplier model:

< λ, u -ḋ > L 2 (Ω O ) = ˆΩO λ • (Π u -Π ḋ) dΩ (4.29)
To partition the energies between the discrete and continuum models, three weight functions α, β and γ are introduced for respectively the internal energy, kinetic energy and the external work of the continuum subdomain Ω C . All these functions verify the following condition (Fig. 4.11):

α, β, γ : Ω C → [0, 1] x →    1 in Ω C \Ω O [0, 1] in Ω C|Ω O 0 otherwise (4.30)
In a complementary manner, the internal energy, the kinetic energy and the external work of the discrete subdomain Ω D are weighted by respectively ᾱ = 1α, β = 1β and γ = 1γ.

Global weak formulation

Based on the previous weak formulations of the isolated discrete and continuum subdomains, the coupling conditions can be introduced to obtain the global weak formulation of the whole problem. By introducing the weight functions in (4.25) and (4.27) and the junction model (4.28), the global weighted weak formulation can be derived:

Find (u, d, θ, λ) ∈ U ad × D ad × O ad × M such that ∀ (δ u, δ ḋ, δ θ, δλ) ∈ U ad,0 × Ḋad,0 × Ȯad,0 × M, given the initial conditions:

ˆΩC β ρ δ u • ü dΩ + ˆΩC α ε(δ u) : A : ε(u) dΩ - ˆ∂Ω T C γ δ u • t g dΓ - ˆΩC γ ρ δ u • f dΩ + n p ∑ i=1 βi m i di • δ ḋi + n p ∑ i=1 βi I i θi • δ θi - n p ∑ i=1 ( γi f ext i + ᾱi f int i ) • δ ḋi - n p ∑ i=1 ( γi c ext i + ᾱi c int i ) • δ θi + δ ˆΩO λ • (Π u -Π ḋ) + l 2 ε(λ) : ε(Π u -Π ḋ) dΩ = 0 (4.31)
Since the last formulation (4.31) is true for any infinitesimal variation of u, ḋ, θ and λ, it can be reformulated as:

Find (u, d, θ, λ) ∈ U ad × D ad × O ad × M such that ∀ (δ u, δ ḋ, δ θ, δλ) ∈ U ad,0 × Ḋad,0 × Ȯad,0 × M,
given the initial conditions:

• DEM side:

n p ∑ i=1 βi m i di • δ ḋi - n p ∑ i=1 ( γi f ext i + ᾱi f int i ) • δ ḋi - ˆΩO λ.δΠ ḋ + l 2 ε(λ) : ε(δΠ ḋ) dΩ = 0 n p ∑ i=1 βi I i θi • δ θi - n p ∑ i=1 ( γi c ext i + ᾱi c int i ) • δ θi = 0 (4.32)
• CNEM side:

ˆΩC β ρ δ u • ü dΩ - ˆ∂Ω T C γ δ u • t g dΓ + ˆΩC α ε(δ u) : A : ε(u) dΩ - ˆΩC γ ρ δ u • f dΩ + ˆΩO λ • δΠ u + l 2 ε(λ) : ε(δΠ u) dΩ = 0 (4.33)
• Interface:

ˆΩO δλ • (Π u -Π ḋ) + l 2 ε(δλ) : ε(Π u -Π ḋ) dΩ = 0 (4.34)

Discretization and spatial integration

Using the CNEM discretization, the continuum subdomain Ω C is represented by a set of scattered nodes (with no connectivity between them). As mentioned above, the CNN interpolation is also introduced in Ω D|Ω O so that continuous field variables would be obtained from the discrete quantities defined at the DEM particle positions in Ω O . Consequently, Ω D|Ω O is assumed to be a continuum towards the interpolation and the DEM particles are also considered as CNEM nodes. The CNN interpolation is only applied in Ω D|Ω O which is assumed to be far from the fine scale effects and not in the whole Ω D .

In this study and contrary to previous works on the discrete-continuum coupling, the dual geometric space of M (mediator space) is taken independent from Ω D|Ω O and Ω C|Ω O . Therefore, a third subdomain proper to the overlapping region is introduced which will be designated hereafter by Ω O (overlapping subdomain). Consequently, the discretization of this subdomain can be made independently of the continuum and discrete ones. This makes easier the application of the DEM-CNEM coupling method for complex problems. Also, this provides more flexibility to the discretization quality of the continuum and discrete subdomains in the overlapping region ( §4.5). The overlapping subdomain is also discretized by a set of scattered nodes. The discretized subdomains of Ω D , Ω C and Ω O are denoted by

Ω h D , Ω h C and Ω h O , respectively. Since Ω D is naturally discretized, Ω h D is the same that Ω D (Ω h D = Ω D ).
After the geometries discretization, the discretized form of the associated (dual) vector spaces D ad , U ad and M can be derived. Their discretized forms are denoted by D ad,h , U ad,h and M h , respectively. According to the configuration of the discretized subdomains in the overlapping zone, four cases can be distinguished (Fig.

4.13).

The first three configurations (Figs. 4.13a, 4.13b and 4.13c) present some regularity in the overlapping zone and are widely studied in the literature. In this work, the fourth configuration which is the most general one including the others, is retained to perform the DEM-CNEM coupling. In this configuration no coincidence conditions are imposed in the overlapping region. Therefore, it is more adapted for complex geometries. In this case, it is sufficient to discretize the different subdomains independently and mount them as indicated in Figure 4.13d. In fact, using this configuration, it is very difficult to prove mathematically the existence and uniqueness of the solutions. Also, and contrary to the three other configurations, no numerical works studying this configuration exist in the literature. In this work, the well-posedness of the global problem will be analyzed numerically in Section 4.5.

Using the CNN interpolation on the different discretized vector spaces, D ad,h , U ad,h and M h , the displacement fields d and u and the Lagrange multipliers λ are approximated by: Concerning the integrals resulting from the application of the coupling conditions, this technique cannot directly be applied. Indeed, these integrals include variables of different natures and defined on different Voronoï diagrams. The issue here is how to choose the background of the integration in the overlapping region. In this work, the Voronoï cells associated with the overlapping subdomain Ω h O are chosen as an integration background in this zone. All of the variables that are not defined in M h are projected on this space using the CNN interpolation. After integration, the discretized equations can be derived:

d h (x) = n D ∑ i=1 φ D i (x) d i (4.35) u h (x) = n C ∑ i=1 φ C i (x) u i (4.36) λ h (x) = n O ∑ i=1 φ O i (x) λ i (4.
a) Ω h C|Ω O = Ω h O = Ω h D|Ω O (b) Ω h C|Ω O = Ω h O ⊂ Ω h D|Ω O (c) Ω h C|Ω O ⊂ Ω h O = Ω h D|Ω O (d) Ω h C|Ω O = Ω h O = Ω h D|Ω O
• DEM side m β d = f int α + f ext γ + f λ I β θ = c int α + c ext γ (4.38)
with:

m β,ij = β(x i ) m i δ ij , m i
and x i are the weighted mass and coordinates of the particle i, respectively.

f int α,i = ᾱ(x i ) f int i :
weighted internal forces applied on particle i by its neighbors.

f ext γ,i = γ(x i ) f ext i :
weighted external forces applied on particle i. I β,ij = β(x i ) δ ij I i , I i is the weighted mass moment of inertia of particle i. c int α,i = ᾱ(x i ) c int i : weighted internal torques applied on particle i by its neighbors.

c ext

γ,i = γ(x i ) c ext i : weighted external torques applied on particle i

f λ = ( c L2 D + l 2 c H1 D ) {λ} = [c D ]
{λ}: coupling forces applied on the par- ticles belonging to the overlapping region. c L2 D and c H1 D are the discrete L2 and H1 coupling matrices.

• CNEM side

M β { ü} = -F int α + F ext γ -F λ (4.39)
with:

M β,ij = β(x i ) M i δ ij , M i is the weighted lumped mass of the node i of coordi- nates x i . F int α = [K α ]
{u} is the vector of the weighted internal forces, [K α ] is the weighted stiffness matrix.

F ext γ
is the vector of the weighted external forces.

F λ = ( C L2 C + l 2 C H1 C ) {λ} = [C C ]

Time integration

The matrix systems (4.38), (4.39) and (4.40) are in the form of time-dependent ordinary differential equations (ODEs). These equations will be solved numerically. Many types of time integration schemes can be found in the literature which can be classified into three headings: explicit, implicit and semi-implicit (or also semiexplicit). The explicit schemes calculate the state of a system at a later time from its state at the current time. The implicit schemes use both current and later states to find the solution at a later time. The semi-implicit schemes combine some features of the explicit and implicit schemes to find the solution.

Since the DEM-CNEM coupling approach, developed in this work, is mainly designed for dynamic simulations, the class of explicit schemes is retained in this work. This class is more adapted for highly dynamic problems such as the Laser Shock Processing. Many explicit schemes can be found in literature, such as the Runge-Kutta, position and velocity Verlet schemes. A comparison between these schemes can be found in Reference [START_REF] Rougier | Numerical comparison of some explicit time integration schemes used in DEM, FEM/DEM and molecular dynamics[END_REF]. According to this reference, the Verlet scheme provides good results. Also, it is easier to implement. Therefore, this scheme is retained to solve the above matrix systems. This scheme gives an O(∆t 3 ) approximation for both velocities and displacements. Thus, the choice of a velocity coupling (control of velocities in the overlapping region) does not affect the coupling approach accuracy compared with the displacement coupling (control of displacements in the overlapping region).

Algorithmic

Now, the different ingredients for a numerical computation are well-designed and defined. Therefore, the associated algorithms which will be translated into computer codes in C++ and python languages can be derived (Alg. 1, 2 and 3). The key idea of these algorithms is: at each time step, ignoring the coupling conditions, predictive accelerations and velocities are determined for the two models (DEM and CNEM); these quantities are then sent to a DEM-CNEM Interface (DC-Interface) code to be corrected taking into account the coupling forces; finally, the corrected quantities will be returned to the discrete and continuum codes to continue computation. This process is repeated until the simulation is achieved.

DEM resolution

• Initialization d n , . d n , {d} n , θ n ,
. θ n and {θ} n : initial conditions or the DC-Interface results.

• Computation of {d} n+1 and {θ} n+1 using the Verlet velocity scheme:

{d} n+1 = {d} n + ∆t . d n + ∆t 2 2 d n {θ} n+1 = {θ} n + ∆t . θ n + ∆t 2 2 θ n (4.41) • Computation of f int α n+1 , f ext γ n+1
, c int α n+1 and c ext γ n+1

• Computation of the predictive linear accelerations d * n+1 (omitting the coupling forces f λ from Equation 4.38).

d * n+1 = m β -1 ( f int α n+1 + f ext γ n+1 ) (4.42)
• Computation of the angular accelerations θ n+1 : u n and {u} n : initial conditions or the DC-Interface re- sults.

θ n+1 = I β -1 ( c int α n+1 + c ext γ n+1 ) (4.
• Computation of {u} n+1 :

{u} n+1 = {u} n + ∆t . u n + ∆t 2 2 { ü} n (4.46)
• Computation of the predictive linear accelerations { ü} * n+1 : (omitting the coupling forces F λ from Equation 4.39).

{ ü} * n+1 = M β -1 (-F int α n+1 + F ext γ n+1 ) (4.
u n+1 = . u * n+1 -∆t 2 M β -1 F λ n+1 . d n+1 = . d * n + ∆t 2 m β -1 f λ n+1 (4.49) F λ n+1 = [C C ] {λ} n+1 f λ n+1 = [c D ] {λ} n+1 (4.50) [C O ] . u n+1 -[c O ] . d n+1 = 0 (4.51)
Introducing Equations (4.49) and (4.50) into Equation (4.51), the interface system of equations can be written as:

[A] {λ} n+1 = {b} n+1 (4.52) where the coupling matrix [A] and {b} n+1 are respectively defined as: • The linear acceleration corrections: { ü} n+1 and d n+1 :

[A] = ∆t 2 ([C O ] M β -1 [C C ] + [c O ] m β -1 [c D ]) (4.53) {b} n+1 = [C O ] . u * n+1 -[c O ] . d * n+1 (4.
{ ü} n+1 = { ü} * n+1 -M β -1 F λ n+1 d n+1 = d * n+1 + m β -1 f λ n+1 (4.
{ u} n+1 { ü} n+1 → CNEM code n ← n + 1 end for

Remarks

• To compute the predictive accelerations, the lumped mass matrices must be invertible. Thus, the weight functions β and β must be strictly positive in Ω O and at the border ∂Ω O . A small ε will be used instead of zero at the nodes assigned to ∂Ω O . The definition of this weight function β given by (4.30) is slightly modified as:

β : Ω C → [0, 1] x →    1 in Ω C \Ω O [ε, 1 -ε] in Ω C|Ω O 0 otherwise (4.56)
where ε is a small strictly positive real number to be chosen.

• To determine the Lagrange multipliers vector at a later time {λ} n+1 , the res- olution of the matrix system (4.52) is performed using the well-known LU decomposition method [START_REF] Meyer | An efficient implementation of LU decomposition in C[END_REF]. In this system, the coupling matrix [A] (4.53) can directly be determined. Since, the mass matrices of both discrete and continuum models are diagonal, it is therefore easy to determine their inverse matrices and to compute [A].

• From an algorithmic point of view, the velocity coupling used in this work is easier than the displacement coupling. Indeed, the displacement coupling requires the computation of the predictive displacements, in addition to the predictive accelerations and velocities, which must be sent to the DC-Interface code for correction.

• No coupling conditions are explicitly applied to correct the angular velocities and angular accelerations of the particles in the overlapping zone. These quantities are corrected implicitly through the correction of the linear velocities and accelerations. Indeed, the internal forces are computed accounting for particle displacements and rotations. These forces are later used to compute the new displacement which are corrected using the coupling condition (4.40).

Implementation

In this work, two existent workbenches have been retained to perform the DEM and CNEM computations which are respectively GranOO (Granular Object Oriented) and NESSY (NEtwork Solver SYstem). GranOO has been developed at the Mechanics and Engineering Institute of Bordeaux (I2M) by J. L. Charles et al. This code provides C++ Object Oriented (OO) libraries to describe and solve dynamic mechanical problems by DEM using explicit temporal integration schemes. This workbench enables customization by adding plug-ins. To implement the DEM-side coupling algorithm in this code, several plug-ins are written in C++ language using the advanced OO architecture.

NESSY has been developed at the Process and Engineering in Mechanics and Materials laboratory (PIMM) by P. Lorong et al. This workbench provides C++ libraries interfaced with Python moduli. The CNEM-side coupling algorithm is therefore implemented in Python language.

The communication between the two codes is ensured by "DC-Interface" developed in this work. DC-Interface which is written in Python language communicates directly with NESSY using Python classes. Whereas, the communication between DC-Interface and GranOO is performed by an Inter Process Communication (IPC) technique. IPC techniques can be divided into mainly four classes: Message Passing, Inter Process Synchronization, Shared Memory, and Remote Procedure Calls. In this work, the Inter Process Synchronization is retained. Using this technique, the synchronization of the different codes in communication is automatically ensured, which is computationally beneficial.

Parametric study of the coupling parameters

Several works have studied mathematically the Arlequin method for both continuumcontinuum coupling [START_REF] Rateau | Méthode Arlequin pour les problèmes mécaniques multi-échelles Applications à des problèmes de jonction et de fissuration de structures élancées[END_REF][START_REF] Dhia | Analyse mathématique de la méthode Arlequin mixte[END_REF][START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF] and continuum-discrete coupling [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF]. The main results concerning the well-posedness of the coupling problem are recalled hereafter. The weight function α, applied to weight the internal energy, must be strictly positive in the overlapping zone. Without this condition the coercivity of the internal energy cannot be verified. Another significant result concerning the coupling junction models is that for the discretized problem, contrary to the H1 coupling (based on the H 1 (Ω O ) scalar product) which yields a well-posed problem, the L2 coupling (based on the L 2 (Ω O ) scalar product) can lead to an ill-conditioned system of equations, especially in the case of highly refined discretization. In the same context, Bauman et al. [START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF] have studied another coupling model, the H1 semi-norm coupling, in which the first term of the H 1 (Ω O ) scalar product (4.28) is removed. This model leads to a well-posed problem, but it does not constrain enough the continuum and discrete displacements in the overlapping zone. Other works [START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF][START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF] have studied numerically the ingredients of the Arlequin method using 1D models. Guidault et al. [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF] have noted that, for the L2 coupling, the weight function α must be continuous at the boundary of the gluing zone ∂Ω O . Indeed, the use of a discontinuous weight function can cause undesirable free conditions at ∂Ω O . Concerning the choice of the mediator space, Ben Dhia et al. [START_REF] Dhia | Further Insights by Theoretical Investigations of the Multiscale Arlequin Method[END_REF][START_REF] Dhia | Analyse mathématique de la méthode Arlequin mixte[END_REF] mentioned that in the case of continuous domains, it is convenient to choose M = H 1 (Ω O ); however, it is very difficult to choose the finite approximation space M h O . To address this difficulty, several works [START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF][START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF] proposed a 1D numerical study of M h O . The different configurations that were studied are presented in Figures 4.13a, 4.13b and 4.13c. The static studies of Guidault et al. [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF] show that: (i) in the case of a fine mediator space M h O (Fig. 4.13c), the response of the structure do not depend on the weight functions and a locking phenomenon takes place, i.e. the fine solution exactly conforms to the coarse solution in the overlapping zone; (ii) in the case of a coarse mediator space (Fig. 4.13b), the weight functions have an influence on the solutions such that the larger the weight function on the fine mesh, the smaller becomes the maximum jump between the two meshes. This work extends these studies to 3D dynamic studies using the general configuration given in Figure 4.13d. It will be demonstrated that some of the results proven in statics using 1D models are not valid in 3D dynamic simulations.

Assuming the general case of the approximated mediator space (Fig. 4.13d), the various coupling parameters studied in this work are:

• The junction parameter l.

• The weight functions α, β and γ.

• The width of the overlapping region L O .

• The discretization of the approximated mediator space M h O . A 3D beam model is used for the dynamic study (Fig. 4.15), of which the length and the diameter are L = 20 mm and D = 2 mm, respectively. This model is divided into two subdomains with an overlapping zone. The left subdomain is modeled by the CNEM approach using 626 nodes (the associated characteristic length is about l Ω C c = 0.47 mm) and fixed at the left end (x = 0). The right subdomain is modeled by the DEM approach using 20 000 spherical particles having r p = 0.05 mm as mean radius. Based on the discretization characteristic lengths of DEM and CNEM subdomains (l Ω D c ≈ r p and l Ω C c ), the cutoff frequencies of the two coupled models can be determined: f CNEM c = 1.9 MHz and f DEM c = 18.2 MHz. To generate a state of high frequency waves in propagation within the beam and investigate the wave reflexion at the DEM-CNEM interface, the free end (x = L) is submitted to a tensile loading with a very steep slope (Fig. ). The material of the beam is the silica, whose the mechanical properties are given in Table 2.1 and recalled in Figure 4.15. The corresponding microscopic properties of the cohesive beam bonds in the DEM approach are given in Table 4.1. For the sake of simplicity, the same weight function is used to weight the different energies (α = β = γ). To control the wave propagation within the model, four check points are placed along this beam (Fig. 4.15) as follows:

• CnemCheckPoint: at the middle of the CNEM subdomain where the controlled quantities are computed using the CNEM nodes in this zone

• OverlapCnemCheckPoint: at the middle of the overlapping zone where the controlled quantities are computed using only the CNEM nodes in this zone.

• OverlapDemCheckPoint: at the middle of the overlapping zone where the controlled quantities are computed using only the DEM particles in this zone.

• DemCheckPoint: at the middle of the DEM subdomain where the controlled quantities are computed using the DEM particles in this zone.

Figure 4.17 presents the reference results obtained by DEM and CNEM separately. Table 4.2 presents the mean displacement of the right end and the first three natural frequencies. It can be seen that the results are in good agreement, and they are also in agreement with the beam theory results. This ensures the equivalence of the two models.

Influence of the junction parameter l

This parameter is mainly employed to compute the coupling matrix [A] (4.53). The influence of this parameter on the coupling results will therefore be investigated by studying its effect on the conditioning of [A] defined by:

Cond([A]) = [A] . [A] -1 (4.57)
which is an important indicator of the results accuracy. A coarse mediator space M h O (coarse Ω h O , i.e. at the same order as Ω h C ) and continuous weight functions with ε = 0.05 (ε is defined in (4.56)) are used in this study. Figure 4.18 shows the conditioning of [A] with respect to l. The [A] conditioning decreases with l and reaches a minimum at a small value (l = l opt ). Beyond this value, the conditioning increases exponentially as l increases. Therefore, the H1 coupling (4.28) for a small value of l Figure 4.17: Free-end displacements obtained using DEM and CNEM separately and the associated spectral analyses (computed from FFT) (width and discretization characteristic length of Ω h O , respectively). l opt practically does not change with the width of the overlapping region, except for very fine mediator space (case of l Ω O c = 32 µm). However, this parameter depends significantly on l Ω O c . In the case of relatively coarse mediator space (case of l Ω O c = 320 µm), l opt is roughly at the same order of magnitude as l Ω O c . In practice, the junction parameter l can be chosen as: l = l Ω O c . This choice generally results in acceptable conditioning.

U mean (mm) f 0 (Hz) f 1 (Hz) f 2 (

Influence of the weight functions α

The influence of the weight function α will be studied in this subsection. The other coupling parameters are chosen as follow: a fine mediator space M h O (fine discretization of Ω h O , i.e. at the same order as the discretization of Ω h D ), L O = 2 mm and α = β = γ.

4.20).

This figure shows that the High Frequency Waves (HFWs) initially captured at the "DemCheckPoint" did not appear in "OverlapDemCheckPoint" or "OverlapC-nemCheckPoint". In other words, the major part of the HFWs are reflected without entering the overlapping region. This explains the deviation in the temporal displacement each time the global wave crosses the overlapping region. Thus, constant weight functions are not a good choice for dynamic simulations. Indeed, the projection mechanism, which occurs in Ω O , cannot dampen the HFWs, and additional filtering is required in this case. In contrast, in the case of static studies, Guidault et al. [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings[END_REF] have shown that constant weight functions can be used with H1 coupling. A large difference between the results is observed. In the first case (α = 0.3), the magnitude of the free-end displacement is greater than the one obtained using α = 0.5. However, it is smaller in the case of α = 0.8. To provide an explanation for these results, the temporal velocities at the check points are presented in Figure 4. [START_REF] Benson | Computational methods in Lagrangian and Eulerian hydrocodes[END_REF].

For α = 0.8, a portion of the global incident wave is positively reflected at the interface between the two models without entering the overlapping zone. Furthermore, only a complementary part is transmitted in the CNEM model. Quantitatively, the transmission and reflection coefficients are t num = 0.44 and r num = 0.56, respectively. By analogy with the wave propagation between media with different acoustic impedances, we define the transmission and reflection coefficients in terms of α and ᾱ as follows: t th = 2 ᾱ α+ ᾱ and r th = α- ᾱ α+ ᾱ . It can be verified that t num and r num are of the same order of magnitude as t th = 0.4 and r th = 0.6, respectively. For α = 0.3, the same reflection mechanism takes place but with a negative coefficient. Indeed, the velocity magnitude of the transmitted wave (measured at "Cnem-CheckPoint") is greater than the velocity magnitude of the incident wave (initially measured at "DemCheckPoint"). Also, the numerical transmission and reflection coefficients, t num = 1.44 and r num = -0.44, are of the same order of magnitude as t th = 1.4 and r th = -0.4. Then, for the case of constant weighting, the use of α = 0.5 is worse with regard to the reflection phenomenon. In this case, a part of the global incident wave will be reflected at the interface between the models (not only the HFWs which are not supported by the coarse scale model). This result proves that the 1D static studies available in the literature cannot be used to perform dynamic coupling. Guidault et al. [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF] have shown that, in statics and using a fine mediator space, the solutions do not depend on the weight functions.

Continuous weight functions

As explained in Remark 2 ( §4.4.5), the weight functions must not vanish at the boundary of the overlapping zone ∂Ω O , and a small value ε must be applied rather than zero at ∂Ω O . Before studying the influence of the continuous weight functions, the influence of ε is investigated. Figure 4.23 presents the free-end displacement using continuous weight functions for different ε. This parameter ε, when less than 0.05, has no practical influence on the results, but a very small ε can lead to instability problems. Indeed, as shown in Table 4.4, the smaller the ε, the greater the conditioning of the coupling matrix [A] becomes. In the remainder of this manuscript, ε = 0.05 will be used each time a continuous weight function is used. Figure 4.24 shows the free-end displacement for the case of a linear continuous weight function (Fig. 4.11). The associated velocities measured at the different check points are presented in Figure 4.25. No high frequency waves (HFWs) are reflected at the interface between the two models. Using continuous weight functions, the HFWs enter the overlapping zone, and then, they are dampened by the projection of the fine space solution (from the DEM model) onto the coarse space solution (from the CNEM model). Figure 4.25b evidences that with a fine mediator space, a small overlapping zone is sufficient to cancel out all of the HFWs. The use of a continuous weight function significantly improves the results (Fig. 4.24a). However, a small deviation from the reference results still persists and becomes greater each time the wave travels back (CNEM-DEM direction). Because of the very fine discretization of the DEM subdomain, the weight of the particles in the overlapping zone decreases smoothly when approaching the CNEM domain.

By examining the first round-trip in Figure 4.24a, it is apparent that practically no deviation from the reference is noted when the wave travels from the DEM subdomain to the CNEM subdomain. Therefore, the incident wave correctly crosses the interface between the pure DEM (Ω D \Ω O ) and the overlapping zone (Ω O ). However, in the CNEM subdomain, a coarse discretization is used. The jump between the weights of the adjacent nodes is relatively large. Thus, the same reflection mechanism, observed using constant weight functions with α = 0.5, occurs when the wave travels back (CNEM-DEM direction). To enhance the results, one can increase the width of the overlapping zone to reduce the slope of the weight functions and so the weight jump between the adjacent nodes. Another solution consists in using continuous differentiable weight functions (Fig. present the results using the two solutions. The wave correctly crosses Ω O without any significant deviation. In conclusion, the continuous differentiable weight functions (Fig. 4.11) are the best-adapted functions to weight the different energies in the overlapping zone. These functions ensure a better regularity (smoothing) in the vicinity of the ∂Ω O , and then the HFWs cross the overlapping zone where they will be dampened, provided that L O and M h O are well-chosen. Otherwise, the same reflection mechanism can take place but for other reasons as will be explained in the next subsections. In the previous subsection, a fixed fine mediator space has been used. In such a case, the fine scale model velocities are practically locked on those of the coarse scale model (Fig. 4.26a). Indeed, the velocity curve at the "OverlapDemCheckPoint" is superimposed on that at the "OverlapCnemCheckPoint". The same locking phenomenon is noted in the literature when the third configuration (Fig. 4.13b) is used [START_REF] Guidault | On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers[END_REF]. Now, to study the influence of M h O on the coupling results, a coarse mediator space is used, i.e. at the same scale as M h C . As shown in Figure 4.26b, the velocities equality in the overlapping region is satisfied only in a weak sense and not in each node of the mediator space. This allows the fine model (DEM model) to correctly act in Ω O . However, in this case a small overlapping zone is insufficient to correctly transmit the global incident wave and cancel the HFWs. A solution to reduce the wave reflection at the interface between the two models is to use a large overlapping zone (large L O ). As a consequence of this solution, a large discrete subdomain has to be used which is costly and cumbersome. However, contrary to previous works on coupling approaches, the discretization of the mediator space M h O can be made independently from the discrete and continuum models in this work. Therefore, an intermediate discretization of M h O , between the discretization of the continuum and discrete models, can be used. As it will be seen in the next subsections, for any choice of L O , there is an optimal discretization of the mediator space that cancels the HFWs in the overlapping region. O , the width of the overlapping region has no significant influence on the results because of the locking phenomenon. In such a case, the DEM particles are strongly constrained in the overlapping region which in turn will represent a rigid wall for the HFWs, regardless of L O . However, in the case of coarse M h O , the discrete particles can somewhat freely act in the overlapping region. Therefore, using a small L O (narrow overlapping region), the interface between the two models will act like a free surface for the HFWs. In this case, a large overlapping region is required to sufficiently constrain the DEM particles in the overlapping region and then to reduce the HFWs reflection. Figure 4.28 shows that, in the case of coarse mediator space, the coupling results are better as L O increases.

Ω O c is l Ω O c = 2 l Ω D c
. Indeed, no significant deviation is pointed in the associated result. However, in- The same conclusions and remarks, as for L O = 2 mm, are pointed. However, in the present case, the optimal value is l

Ω O c = 9 l Ω D c
and the deviation mechanisms are much less sensitive to l Ω O c . In conclusion, the optimal value of l Ω O c depends also on L O , but the choice of this parameter is less sensitive when large overlapping region is used.

How to choose the coupling parameters in practice?

In a general case there is not an obvious method to determine, in a single way, the various coupling parameters to avoid wave reflection. This subsection gives several recommendations and trends to choose correctly these parameters. To minimize the conditioning of the coupling matrix [A], the H1 coupling with l = l opt is recommended. l opt can be chosen as the discretization characteristic length of the overlapping subdomain (l opt = l Ω O c ). The weight functions must have at least C 0 continuity in the overlapping region, but more regular functions ( having C 1 continuity or more) are recommended. Indeed, with constant weight functions, the high frequency waves (HFWs) are reflected without entering the overlapping region and cannot be dampened by the projection mechanism which occurs in this region. The choice of L O and M h O is more challenging. Indeed, these parameters depend on each other, and on the discretization of the discrete and continuum models. In practice, the choice of these parameters can be made as follows: the discrete and continuum models are discretized according to the physical effects that have to be captured; then, the width of the overlapping region can be fixed as narrow as possible, containing sufficient number of CNEM nodes and DEM particles; finally, M h O can be chosen such that l

Ω D c < l Ω O c < l Ω C c . Using an arbitrary choice of l Ω O c banded between l Ω D c and l Ω C
c , a dynamic simulation can be performed to control the wave propagation between the models; if an inward deviation is found when the wave crossing the interface, l Ω O c must be increased; if an outward deviation is found, l Ω O c must be decreased; Otherwise, the optimal value of l Ω O c is luckily used.

Validation

The previous parametric study using 3D beam submitted to a dynamic tensile loading has highlighted several recommendations to perform a correct coupling. In this Bending U y mean (mm) 5.85e -1 5.86e -1 6.27e -1 6.13e -1 f 0 (Hz)

1 606 1 594 1 557 1 595 Torsion θ x mean (mm) 2.05e -3 2.34e -3 - 2.32e -3 f 0 (Hz) 9 382 9 252 - 9 106
Table 4.5: Comparison of results; L O = 10 mm, coarse mediator space , continuous differentiable weight functions, ε = 0.05 section, the results of this study are used to validate the coupling between CNEM and DEM. Contrary to the tension case, in bending and torsion, the strains in the cross sections are significant. To account for these effects, new geometric characteristics of the 3D model are used: L = 100 mm and D = 20 mm (L/D = 5). The DEM method is applied only for the portion located 20 mm from the right end (the section located at x = L) and the remainder of the domain is modeled using the CNEM method (Fig. 4.32). The different Arlequin parameters are chosen as follows:

L O = 10 mm, continuous differentiable weight functions, ε = 0.05, l = l Ω O c
and coarse mediator space. Figures 4.33 and 4.34 present the temporal free-end displacements with respect to the xaxis and yaxis using respectively tensile and bending loading. The deviation from the reference, as observed in the previous simulations when the wave crosses Ω O , disappeared in the present results. Figure 4.35 presents the temporal free-end rotation about the xaxis using the torsional loading. The coupling results are in good agreement with the DEM results (Tab. 4.5). Finally, the coupling method has been tested using an initial velocity loading (Test 4 in Figure 4.32). Figure 4.36 presents the associated free-end displacement with respect to xaxis. The coupling result is comparable to the reference one. The comparison between the coupling results and the results obtained using DEM and CNEM separately (Figs.

Conclusion

In this chapter, a 3D discrete-continuum coupling method adapted for dynamic simulations has been developed. This approach couples the two methods that have been retained in the previous chapter: the Discrete Element Method (DEM) and the Constrained Natural Element Method (CNEM).

After recalling the main features and specificities of these two methods, the coupling approach has been detailed. This approach is based on the Arlequin technique which gives a flexible framework to couple dissimilar methods. The most general configuration of the superposed subdomains in the overlapping zone Ω O (Fig. 4.13d) was retained to perform the DEM-CNEM coupling. In this configuration no coincidence conditions are imposed in Ω O . Therefore, the different subdomains can be discretized independently, which makes easier the application of this coupling approach on complex geometries. As shown, the well posedness of the global problem was verified numerically.

Several parameters are involved in this coupling approach: the junction parameter l, the weight functions α, β and γ, the width of the overlapping region L O and the approximated mediator space M h O . To simplify the setting of these parameters in practice, a parametric study of these parameters were performed using 3D beam model subjected to dynamic tensile loading. For the case of a small junction parameter l, the H1 coupling is more accurate than the L2 coupling. Indeed, it decreases the conditioning of the coupling matrix [A]. However, beyond a small value (l opt ), it can lead to instability problems. Therefore, it is important to carefully choose this parameter. In practice, l can be chosen as the discretization characteristic length of the overlapping subdomain (l opt = l Ω O c ). This choice generally results in acceptable conditioning. The weight functions must have at least C 0 continuity in the overlapping region, but more regular functions (having C 1 continuity or more) are recommended. Indeed, with constant weight functions, the high frequency waves (HFWs) are reflected without entering the overlapping region and cannot be dampened by the projection mechanism which occurs in this region. The choice of L O and M h O is more challenging. These two parameters depend on each other, and on the discretization of the discrete and continuum subdomains. In practice, L O can first be fixed such that the overlapping region contains sufficient number of DEM elements and CNEM nodes, and then preliminary dynamic simulations using simple dynamic loading and different M h O can be performed to choose the appropriate one. Since the CNEM approach is very close to the FEM method and has practically all its advantages, the performances of the developed coupling approach should not be affected, compared with a DEM-FEM coupling approach. In addition, since the CNEM is a meshfree method, the DEM-CNEM coupling approach is better suited to simulate complex problems. This chapter has been the subject of a first paper published in Computer Methods in Applied Mechanics and Engineering [START_REF] Jebahi | 3D coupling approach between discrete and continuum models for dynamic simulations (DEM-CNEM)[END_REF].

After developing the numerical tool which will be used to simulate the LSP processing on silica glass, it is necessary to develop a silica glass mechanical behavior model adapted for this approach to correctly predict its response to laser-induced pressures. This is the subject of the next chapter.

Chapter 5

Modeling of the silica glass mechanical behavior by the Discrete Element Method

"A theory is something nobody believes, except the person who made it. An experiment is something everybody believes, except the person who made it."

Albert Einstein

Abstract

This chapter is devoted to the numerical modeling of the silica glass mechanical behavior using the DEM-CNEM coupling approach. Since the continuum approach (CNEM) will be applied far from the high loaded regions, only elastic behavior is considered in the continuum subdomains. However, a complex model is proposed to take into account the different phenomena of the silica glass mechanical behavior in the discrete subdomains. Hereafter, we will detail how the silica glass mechanical behavior is modeled using the Discrete Element Method.

Introduction

As seen in Chapter 2, silica glass has a complex mechanical behavior. Furthermore, some of its mechanical properties such as the Young's modulus, Poisson's ratio and density show anomalous dependence on the pressure. Under hydrostatic pressure, this material behaves in a perfectly elastic manner up to the densification pressure P d . The elastic behavior of silica glass can be divided into two parts: nonlinear elastic behavior until a limit noted by P nle , and linear elastic behavior beyond. After P d , silica glass begins to exhibit signs of permanent deformation, which is termed densification, until the saturation pressure P s . Beyond P s , this glass saturates and returns purely elastic. In the literature, different values are given for the bounding parameters of the silica glass mechanical behavior as mentioned in Table 2.2. In this work, the Ji's results [90] are used to perform the static studies. These results are the most recent (more accurate equipments were used to obtain them). Also, Ji [90] has conduced a complete experimental study from the investigation of silica glass behavior to the Vickers micro-indentation of this material using samples manufactured in the same conditions which can reduce the results scattering. As for the dynamic studies, the Marsh's [START_REF] Marsh | LASL Shock Hugoniot Data[END_REF] results are used, which are in agreement with those obtained by Yann [START_REF] Michel | Phénomène d'impact à haute vitesse sur cibles minces fragiles : application au projet de laser mégajoule et à la problématique des débris spatiaux[END_REF], Wackerle [START_REF] Wackerle | Shock wave compression of quartz[END_REF], Sugiura [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF] and Lalle [START_REF] Lalle | Utilisation de l'eau, de la silice fondue et de deux verres de la société Shott comme fenêtre[END_REF]. In this chapter, we shall try to model the different phenomena characterizing the silica glass mechanical behavior. This part concerns only the discrete approach (DEM) which is deemed to be applied in the high loaded regions. However, only an elastic constitutive law will be applied in the continuum subdomains (CNEM regions).

A densification model based on the normal stress in the cohesive beams has been developed. This model involves several microscopic parameters which must be determined by calibration. To facilitate the calibration process, the influence of each microscopic parameter has been studied. The validation of this model has first been performed in statics by simulation of static Vickers micro-indentation, and then in dynamics by simulation of high velocity impacts of silica glass plates. Although this model gives relatively good results, the large number of the involved microscopic parameters makes tedious the calibration process. To overcome this limitation, another model based on the virial stress at the discrete elements [START_REF] Mclellan | Virial theorem generalized[END_REF][START_REF] Zhou | A new look at the atomic level virial stress: on continuum-molecular system equivalence[END_REF] has been developed. This model involves only one parameter to be calibrated. The validation of this model has been performed in dynamics (the scope of this work) by simulation of high velocity impacts of silica glass plates. The associated results are very close to those obtained using the first model. Then, it is retained to simulate the LSP processing on silica glass (Chapter 6).

To model fracture by discrete approaches, the models most commonly used are based on the computation of equivalent stresses or strains in the cohesive bonds. Such models give relatively good results at the macroscopic scale. However, they do not reproduce correctly the cracking pattern at the microscopic scale. To solve this problem, the notion of virial stress has also been used to develop a new fracture model. Using this model, relatively good results have been obtained at both macroscopic and microscopic scales.

Modeling hypotheses

Almost all the works on densification of silica glass agree that this phenomenon can take place under hydrostatic pressure. However, the shearing effects on this phenomenon present a central issue for several studies and need more investigations. The works on this subject which can be found in the literature are rather controversial. Indeed, several researchers [START_REF] Cohen | Effects of Ultra high Pressures on Glass[END_REF][START_REF] Cohen | Reply to "Comments on'Effects of Ultrahigh Pressures on Glass[END_REF][START_REF] Dachille | High Pressure Studies of the Silica Isotypes[END_REF] state that shear has only a kinetic effect and does not change the densification level of silica glass (it only changes the rate at which the silica glass responds to pressure), while others [START_REF] Rouxel | Indentation deformation mechanism in glass: Densification versus shear flow[END_REF][START_REF] Christiansen | Irreversible Compressibility of Silica Glass as a Means of Determining the Distribution of Force in Highpressure Cells[END_REF][START_REF] Mackenzie | High-Pressure Effects on Oxide Glasses: I, Densification in Rigid State[END_REF][START_REF] Weir | Comments on "Effects of Ultrahigh Pressures on Glass[END_REF] hold the opposite view. In this work, it is assumed that the shearing effects on the densification are insignificant and they are dominated by the hydrostatic pressure effects. Therefore, densification of silica glass will be supposed to occur only under high hydrostatic pressures. There are many physical and mechanical reasons that support this assumption:

• According to the characteristics of materials, two main mechanisms can explain a permanent deformation due to shear stresses: (i) in the case of crystalline solids, plasticity is due to dislocations slip, defined as crystallographic defects within the crystalline structure. Silica glass as a typical amorphous material lacks the long-range order characteristic of a crystal. Therefore, it is very difficult for dislocations to develop and move within its irregular structure [START_REF] Kanel | Some new data on deformation and fracture of solids under shockwave loading[END_REF]. (ii) In the case of brittle materials, such as rock and concrete, permanent deformation is caused predominantly by slip at micro-cracks. This generally weakens the considered material and degrades its mechanical properties [START_REF] Lubliner | Plasticity theory[END_REF]. This is not the case for silica glass. In fact, permanent deformation in silica glass enhances the mechanical properties [START_REF] Keryvin | Contribution à l'étude des mécanismes de déformation et de fissuration des verres[END_REF]90].

• Because silica glass is an isotropic material, permanent deformation due to shear stress is volume-conservative. Therefore, even if a shear stress can cause silica glass plasticity, this does not affect the densification level (permanent volume change at a given pressure). For an infinitesimal strain, the volume change can be expressed as:

dV V 0 = J -1 ∼ = trace(ε) (5.1)
where dV is the volume change, V 0 is the initial volume, J is the determinant of the Jacobian matrix [START_REF] Wu | Continuum Mechanics and Plasticity[END_REF] and ε is the strain tensor. If only shear stresses are applied, trace(ε) = 0, and therefore dV = 0.

Based on the previous assumption (the shearing effects on the densification are negligible compared with those of hydrostatic pressure), two densification models adapted for discrete methods have been developed in this work. As will be seen, comparison between experimental and numerical densification results using these models has yielded important conclusions on the validity of this assumption.

Concerning the fracture of silica glass, it will be supposed that completely brittle fracture occurs in this material. As reported in several previous works [START_REF] Griffith | The Phenomena of Rupture and Flow in Solids[END_REF][START_REF] Wiederhorn | Fracture of silicate glasses: ductile or brittle?[END_REF], this fracture type is thought to be initiated under tensile stress in mode I.

Mechanical behavior modeling based on normal stress in the cohesive beams

The DEM variant used in this work models a continuum by a set of spherical rigid particles linked by cohesive beam bonds which model the mechanical behavior of the studied material. Therefore, one approach to model the different stages of the silica glass behavior is to introduce a complex mechanical behavior in these beams. The present section details how this aim is put into practice. To ensure the bridging between the macroscopic and microscopic (at the cohesive beams level) behaviors, a simplifying assumption is used. This assumption consists in considering that the hydrostatic pressure in a discrete element i is only determined by summation of the normal stresses in the beams connected to this element (5.2).

P i hyd = n v ∑ j=1 σ n, ij µ (5.2)
where n v is the total number of the cohesive beams connected to the discrete element i and σ n, ij µ is the normal stress in the beam connecting the discrete element i and j. Note the normal stress in a cohesive beam can easily be determined, based on the Euler-Bernoulli beam theory.

Modeling of nonlinear elasticity

For relatively low pressures (under approximately P nle = 3 GPa (Fig. 2.15)), the silica glass shows a nonlinear elastic behavior which translates into pressure dependence of the bulk and shear moduli noted by K and G, respectively. Based on Figure 2.15 which gives the pressure dependence curves of these moduli, the relationships between the pressure and macroscopic elastic moduli (E and ν) can be determined, using (5.3). Figure 5.1 gives the variations of E and ν with pressure.

E = 9 K G 3 K + G and ν = 3 K -2 G 2(3 K + G) (5.3)
To take into account the macroscopic mechanical properties dependence with pressure in the numerical model of silica glass, it is necessary to find the variations in the microscopic cohesive beam properties which imply such macroscopic variations. In other words, it is necessary to find relationships between the macroscopic and microscopic mechanical properties. Based on the work of André et al. [3,4], the following conclusions can be drawn:

• The microscopic Poisson's ratio of the cohesive beams ν µ has no influence on the macroscopic elastic behavior.

• The macroscopic Poisson's ratio depends only on rµ . The function relating these two parameters may be approximated by a second-order polynomial.

• The macroscopic Young's modulus E is a function of both rµ and E µ . For a given value of rµ , E depends linearly on E µ , whereas it is a second-order

Modeling of densification

To model the densification of silica glass a "spring-slider" system is connected in series with each cohesive beam linking two particles (Fig. 5.2). When the compression stress in the associated beam σ c µ reaches the microscopic densification pressure σ d µ , the slider activates and the densification mechanism takes place. The densification level beyond σ d µ is controlled by the microscopic tangential modulus E t µ , which is adjustable (Fig. 5.3). Finally, to model the saturation stage, the slider stops slipping above a certain limit controlled by the maximum microscopic permanent deformation of cohesive bonds ε p max µ . For the sake of simplicity, the variations of the silica glass macroscopic properties (Young's modulus E and Poisson's ratio ν) with hydrostatic pressure (Fig. 2.16) are approximated in this work by linear piecewise functions (Fig. 5.4). The mechanical macroscopic properties remain constant until the densification pressure P d is reached, then, increase linearly between P d and P s (the saturation pressure) and remain constant above P s . To introduce the variation of the macroscopic properties with pressure in the densification model, the same technique explained in Section 5.3.1 is used. Using Equations (5.4) and (5.5), the variations of the microscopic properties of the cohesive bonds with the normal stress can be determined and then taken into account in the densification model. Figure 5.5 presents the cohesive bond behavior taking into account the variation of the microscopic mechanical properties with densification. The relaxation slope (current microscopic Young's modulus) increases with the densification state. The present densification model involves three . These parameters have to be determined by calibration tests.

Static calibration and validation

Static calibration In the case of static studies, the nonlinear elastic behavior is not of major importance. Such studies generally focus on the initial and final state of the studied material, e.g. permanent deformation. Therefore, the calibration of the associated microscopic parameter σ nle µ will be performed in the next subsection. The present subsection will only detail the calibration process of the densification parameters. Figure 5.6 shows the calibration discrete domain which is a sphere of 100 mm diameter subjected to hydrostatic compression.

Three microscopic parameters are involved in the densification model: the microscopic densification pressure σ d µ , the microscopic tangential modulus E t µ and the maximum permanent deformation of the beam bonds ε p max µ . These parameters will be determined by calibration to obtain the macroscopic behavior. In the case of silica glass, this behavior is characterized in statics by three macroscopic parameters ( §2.2.1): the densification pressure P d , the saturation pressure P s and the volume change ratio ∆V V 0 (or density change ratio ∆ρ ρ 0 ). To simplify the calibration process, the influence of each microscopic parameter on the macroscopic ones is studied.

One major problem of the discrete methods is the dependence of results on the number of the discrete elements (n p ) used to discretize the numerical sample. Therefore, the influence of this parameter on the densification response is first studied.

Macroscopic validation

This paragraph attempts to validate, at the macroscopic scale, the different microscopic densification parameters, obtained by calibration (Tab. 5.1). The same geometric model as for the calibration study is used in this paragraph. The various microscopic parameters are set as indicated in Table 5.1. This model is subjected to various hydrostatic pressures. The associated permanent deformations are measured and compared to the Ji's [90] experimental results obtained using thick samples (Fig. A hemispherical DEM domain of 12 µm radius is used for the Vickers microindentation studies (Fig. 5.12). This domain, whose the spherical part (dome) is fixed, is discretized using n p = 100 000 discrete elements (i.e., the average radius of the discrete elements is around 0.2 µm) to capture the microscopic mechanisms of the densification. It is indented with a square-based pyramidal indenter with included face angles of 136 • (Vickers indenter) and 2 µm in height (Fig. 5.12). This indenter undergoes a moderate displacement in the z direction until the reacting force reaches the desired indentation force. Table 5.2 presents a comparison between the experimental [90,[START_REF] Yoshida | Quantitative evaluation of indentation-induced densification in glass[END_REF] and the present numerical indentation results for low indentation forces (F indent = 0.1 N and F indent = 0.5 N).

For F indent = 0.1 N, the DEM simulation yielded a relatively good estimate of the volume of the indentation print V -and the volume of the piled-up material Figure 5.12: DEM indentation model around the indenter V + . The densified volume V d which is a function of V -and V + compares favorably with the experimental results [90,[START_REF] Yoshida | Quantitative evaluation of indentation-induced densification in glass[END_REF] (in comparison to previous numerical studies [90]). Furthermore, a relatively good numerical result was obtained for the diagonal length of the indentation print (L) which is used to compute the hardness. However, the indentation depth D was low compared with the experimental results [90,[START_REF] Yoshida | Quantitative evaluation of indentation-induced densification in glass[END_REF]. This can be explained by the fact that in the experimental tests, micro-cracks can initiate under the indenter tip (zone of high stress concentration), causing measurement errors in D. The experimental value of D can include not only the depth of the indentation but also the depth induced by the cracks.

For F indent = 0.5 N, the densification state becomes more significant and a far greater densification print is obtained (Fig. 5.14). Compared with the experimental results (Tab. 5.2), the DEM numerical simulation provided a relatively good estimate of the volume of the piled-up material around the indenter V + and the diagonal length L of the indentation print. Concerning D, the analysis of the F indent = 0.1 N case remains valid. However, the gap between the numerical and experimental values of V -is more significant than for F indent = 0.1 N. This can be explained by the side-effects. In the case of F indent = 0.5 N, the densified zone is relatively important compared with the global geometric model. Consequently, the effects of the overstated fixation conditions applied to the dome (the spherical part of the geometric model) become significant. This can increase the region where the densification pressure P d can be reached, hence the overstatement of V -.

One way to reduce the side effects is to use a larger DEM domain. However, such a solution can amplify considerably the computational time, because a very fine discretization must be applied to capture the microscopic densification effects. To circumvent the computation time problem, the DEM-CNEM coupling model is used. The DEM is applied in the indentation region where the densification pres-Figure 5.13: DEM-CNEM coupling indentation model sure (P d ) can be reached, and the CNEM is applied in the remainder of the modeled domain. Figures 5.13 presents the coupling model which is a box of 30 µm sides length. The DEM subdomain is the same as in the previous study (DEM computation). Beyond the indentation region, the CNEM is adopted using 10 000 nodes. Based on the studies of Chapter 4, a 2 µm-thick overlapping zone is used between these two models to ensure correct wave propagation between them. Because the CNEM domain is far from the indentation region, only an elastic behavior is applied in this region.

As shown in Table 5.2, in the case of F indent = 0.1 N, the DEM and DEM-CNEM coupling results are practically the same. However, in the case of F indent = 0.5 N where the side-effects become significant, the DEM-CNEM coupling results are slightly better than those obtained using the discrete model (DEM) only, although they are relatively overestimated compared to the experimental results. Furthermore, the DEM-CNEM coupling improves the shape of the densified zone as shown in Figure 5.15 which presents a cut-away view of the indented model in the case of F indent = 0.5 N. It can be observed that the densified zone obtained from the coupling method (Fig. 5.15c) is closer to Ji's experimental results (Fig. 5.15a) [90] than the den- ). Indeed, it is less deep and broader in the former (coupling result) than in the latter (DEM result). In conclusion, in the case of low indentation forces, the densification model developed in this work gives an acceptable quantitative estimate of silica glass densification at the microscopic scale.

Dynamic calibration and validation

In the previous subsection, the developed densification model is studied in statics to better analyze its behavior at both the macroscopic and microscopic scales and avoid the dynamic effects. The present subsection aims to study this model in dynamics which is the scope of this work.

Dynamic calibration

To calibrate the silica glass numerical model in dynamics, a plate of 10 × 10 mm 2 square base and 3.5 mm thickness is used. The top surface of the plate is subjected to high velocity loadings up to 4000 ms -1 (Fig. 5.16). This geometric model is discretized using n p = 50 000 discrete elements. In order to ensure a state of uniaxial strain, the particles belonging to the lateral boundaries are not allowed to leave their planes.

To model the nonlinear elastic behavior of silica glass in the region of low pressures, only one microscopic parameter is used: the normal stress threshold σ nle µ from which the silica glass behavior returns linear elastic. Therefore, the calibration of this parameter is fairly easy. A series of simulations using different values of σ nle µ has been performed until obtaining the expected macroscopic P nle . Table 5.3 gives the calibrated value of σ nle µ . As for the densification parameters, their calibration must be performed with respect to the Hugoniot curve of silica glass which characterizes the dynamic behavior of this material (Fig. It should be noted that the Hugoniot curve is obtained using impacts of silica glass plates having large lateral dimensions compared to their thicknesses. Therefore, the state of uniaxial stress V -and V + are respectively the volumes of the indentation print and of the piled-up material around the indent, V d is the densified volume, D is the indentation depth and L is the diagonal length of the indentation print.

F indent V -(µm 3 ) V + (µm 3 ) V d (µm 3 ) D (
Figure 5.16: Geometric model for the dynamic calibration of the microscopic nonlinear elastic and densification parameters in the plates center (measuring region) is reduced to a state of uniaxial strain, and the longitudinal stress can be viewed as pressure. Figure 5.17a presents the longitudinal stress dependence of the Young's modulus. This figure can also be regarded as pressure dependence of Young's modulus. For the sake of simplicity, the pressure dependence of the Young's modulus is approximated by linear piecewise function (Fig. 5.17b). E linearly increases between P d et P s and remains constant everywhere else. Table 5.3 gives two proposals of calibrated densification parameters. The first proposal correctly fits the elastic and densification domains of silica glass; however, it overestimates the macroscopic density change ∆ρ ρ 0 (Fig. 5.18). The second proposal corrects this problem, but it loses accuracy at the end of the densification domain as can be seen in Figure 5.18. Since this work aims to simulate the LSP processing on silica using high-power lasers, of which the induced pressures can far exceed P s , the second proposal is retained for this model.

Validation: plates impact To validate this calibrated numerical model of silica glass in dynamics, the high velocity impact tests of silica glass plates, performed by Sugiura et al. [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF], have been reproduced numerically in this work. In these tests, the flayer and target plates are made of the same material which is silica glass (symmetrical impacts). The final particle velocity is expected to be exactly half of the impact velocity, assuming that there is no diminution in particle velocity at any point. Therefore, to reduce the computation time, only a small portion located at the center of the target is considered in this simulation (Fig. 5.19). The validation geometric model is the same as that used for the dynamic calibration. The top surface (impact surface) is subjected to half of the impact velocities (Fig.

5.19).

To reduce the velocity jump at the first time step (for numerical stability requirements), the velocity step functions are changed as shown in Figure 5.20. In order to compare the numerical results with the experimental ones, two measuring points are introduced in the geometric model, which are located at 1 mm and 3 mm distances from the top surface, respectively. Figure 5.21 presents a comparison between the experimental [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF] and numerical velocity profiles corresponding to impact velocity of V Impact = 3310 ms -1 . It should be noted that this the experimental curves show the particle velocity profiles in the measuring points as recorded on the oscilloscopes, and therefore they are given in terms of voltage versus time. Qualitatively, the numerical results are in good agreement with the experimental ones. Figure 5.22 shows a comparison between the numerical and experimental results obtained using impact velocities of 880 ms -1 (elastic region), 2580 ms -1 (densification region: around the HEL) and 3310 ms -1 (densification region: far from the HEL). Of particular observation is the shock fronts (first and second shock fronts) which are not very steep as in the experimental results. This can be explained by the discretization fineness (number of discrete elements used to discretize the probem domain). Indeed, the shock front (even physically) is not an exact discontinuity but it spreads out over several layers of discrete elements. Therefore, a very fine discretization is required to obtain steep shock fronts. In the elastic region (V Impact = 880 ms -1 ), a relatively good estimate of the maximum particle velocity V max p is obtained. However, in the densification region, slight difference between the experimental and numerical values of V max p is observed. The numerical model underestimates V max p . The discrepancy between the numerical and experimental values of V max p increases as the impact velocity increases. This can be explained by the fact that the slope of the impact velocity (V Impact ) steepens up as the latter (V Impact ) increases (Fig. 5.20b). Therefore, additional very high frequency waves (HFWs) exceeding the cutoff frequency of the DEM discretization can be involved. These HFWs which are not supported by the DEM discretization will numerically be filtered, resulting in V max p decrease. Another reason which can explain the V max p understatement is that the equation of state (Hugoniot) is satisfied only at macroscopic scale and not at each discrete element. Therefore, some errors can be introduced in the numerical estimate of V max p , which is computed as the average velocity of the discrete elements belonging to the measuring region. Table 5.4 gives a comparison between the numerical and experimental waves velocities (elastic and plastic waves). The numerical model provides a relatively good estimate of the elastic wave velocity. However, relatively poor estimate is given for the plastic wave velocity. Indeed, fairly similar values of this parameter are given for V Impact = 2580 ms -1 and V Impact = 3310 ms -1 . Besides, relatively great difference is pointed between the numerical and experimental values. The numerical model overestimates this parameter in the case of V Impact = 2580 ms -1 and underestimates it in the case of V Impact = 3310 ms -1 (far from the HEL region). It should be remembered that the velocity of a disturbance in matter is proportionate to ∂P ∂ρ s . Therefore, the similarity between the numerical plastic wave velocities corresponding to V Impact = 2580 ms -1 and V Impact = 3310 ms -1 can be explained by the fact that the calibrated Hugoniot curve (Fig. 5.18) has practically a constant slope in the region of moderate pressures (up to 20 GPa which corresponds to (a) Experimental result from [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF] (b) Numerical result Figure 5.21: Qualitative comparison between numerical and experimental [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF] particle velocity profiles in the region of permanent densification (Impact velocity V Impact = 3310ms -1 ); curves (A) and (B) correspond to the front and rear measuring points, respectively.

V Impact ( ms -1 ) V Impact = 4000 ms -1 , according the silica glass shock polar (Fig. 2.21)). To explain the difference between the experimental and numerical plastic wave velocities, the calibrated Hugoniot curve and that obtained by Sugiura et al. [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF] (whose experimental results are used here) are isolated and zoomed in on the region of densification (Fig.

1

5.23).

Close to the HEL, the slope of the numerical Hugoniot curve is greater than that of Sugiura et al. [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF] (hence, the overstatement of the plastic wave velocity), whereas the opposite trend is shown far from the HEL (hence, the understatement of the plastic wave velocity).

Discussion

In this section, a silica glass numerical model based on the normal stress in the beams has been developed. In addition to the two elastic parameters (E µ and rµ ), this model involves four other microscopic parameters which are the critical nonlinear normal stress σ nle µ , the microscopic densification pressure σ d µ , the microscopic tangential modulus E t µ and the maximum permanent deformation of the beam bonds ε p max µ . These parameters have been determined by calibration. The static validation of this model has been performed at both macroscopic and microscopic scales. Rel-atively good results have been obtained compared with experimental ones [90,[START_REF] Yoshida | Quantitative evaluation of indentation-induced densification in glass[END_REF]. Also, it is validated in dynamics using symmetrical impacts of silica glass plates. In the region of elastic behavior, the numerical results compare favorably with the experimental ones [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF].

However, relatively poor results are obtained in the region of densification. Besides, this model presents some drawbacks. Many microscopic parameters are involved in this model, which makes the calibration process tedious and very timeconsuming. Also, since only the beams in which σ d µ is reached experience the variation of the microscopic properties with normal stress, a state of anisotropy can be established in the material being simulated, depending on the loading type. For example, in the case of uniaxial loading, only the beams in the direction of loading change their properties, whereas the beams perpendicular to this direction keep the initial properties. This problem has not been reported in the case of plates impacts simulations, because such simulations have focused on the 1D shock wave propagation. However, this problem can become annoying in the case of 3D simulations such as the simulation of LSP processing.

To overcome the limitations of this model, another silica glass numerical model has been developed as will be seen in next section. Contrary to the previous beambased model, the new model ensures that the equation of state is satisfied at each discrete element.

Mechanical behavior modeling based on the virial stress in the discrete elements

In continuum approaches, several techniques have been proposed and successfully applied in the literature to model complex mechanical behaviors of materials. Therefore, it can be beneficial to apply these techniques in discrete approaches. To achieve this aim, it is important to know relationships between certain microscopic (of discrete elements) and macroscopic (of continua) quantities, allowing discrete-continuum mechanics bridging. Pioneering work to bridge the discrete and continuum mechanics has been done in this field by Born and Huang [START_REF] Born | Dynamical Theory of Crystal Lattices[END_REF]. They have used an elastic energy approach to evaluate the stress in lattices by means of the Cauchy-Born hypothesis for homogeneous deformation. More recently, another formalism has been developed to measure stress in Molecular Dynamic (MD), the virial stress. This formalism is a generalization of the virial theorem of Clausius (1870) for gas pressure. It has become very popular and widely used in different discrete methods. Therefore, it is retained in this work to bridge the discrete and continuum mechanics.

Virial stress

Principle The virial stress includes two parts [START_REF] Mclellan | Virial theorem generalized[END_REF]: the first part depends on the mass and the velocity (or, in some versions, the fluctuation of the velocity) of particles (discrete elements); the second part depends on inter-particle forces and particle positions, providing a continuum measure for the internal mechanical interactions between particles. In this definition, the average virial stress over a volume Ω around a particle i is given by (5.6).

Π = 1 Ω (-m i ui ⊗ ui + 1 2 ∑ j =i l ij ⊗ f ij ) (5.6)
where m i and ui are respectively the mass and the velocity of the particle i, l ij = l jl i , l i is the position vector of particle i, f ij is the inter-particle force acting on particle i exerted by particle j and ⊗ denotes the tensor product. The virial stress (5.6) is widely used to compute an equivalent of Cauchy stress in discrete systems. Recently, Zhou [START_REF] Zhou | A new look at the atomic level virial stress: on continuum-molecular system equivalence[END_REF] has demonstrated that this quantity is not a measure for mechanical forces between material points and cannot be regarded as a measure of mechanical stress. The same author [START_REF] Zhou | A new look at the atomic level virial stress: on continuum-molecular system equivalence[END_REF] has proposed another formulation to compute the average mechanical stress, including only the second part of (5.6). For a region Ω around a particle i, the average stress is:

σ = 1 2Ω ∑ i∈Ω ∑ j =i l ij ⊗ f ij (5.7)
The expressions (5.6) and (5.7) are developed for Molecular Dynamics (MD), where the inter-particle forces are a function of the Lennard-Jones potential Φ(l ij ) and the distance between the particles l ij = l ij .

f ij = ∂Φ(l ij ) ∂l ij l ij l ij (5 .8)
In this instance, the two expressions of stress (5.6) and (5.7) lead to symmetric tensors. However, this cannot be generalized to all other discrete models. The DEM variant used in this work, where the particles are linked by cohesive beams, does not ensure the symmetry condition of the resulting tensors. A slight modification of (5.7) is then proposed to restore the symmetry condition of the stress tensor as follows:

σ = 1 2Ω ∑ i∈Ω ∑ j =i 1 2 (l ij ⊗ f ij + f ij ⊗ l ij ) (5.9)
Using (5.9), an equivalent Cauchy stress tensor for each discrete element i can be determined as follows (5.10):

σi = 1 2Ω i n v ∑ j=1 1 2 l ij ⊗ f ij + f ij ⊗ l ij (5.10)
where n v is the number of the discrete elements connected to i, Ω = Ω i is a volume associated to i (generally chosen as the volume of the discrete element i).

Discrete-continuum equivalence Since the modeling of silica glass mechanical behavior is resting on the assumption that the shearing effects are negligible, the discrete-continuum equivalence will be studied in terms of pressures. First, a con-stant Ω i is considered, equal to the volume of the associated discrete element i multiplied by the local volume fraction. Using the model of Figure 5.16 with only elastic behavior, the variation of the macroscopic pressure (P hyd ) with respect to the virial pressure (P virial ) in the measuring region is plotted (Fig. 5.24). The choice of constant Ω i leads to nonlinear dependence between P hyd and P virial , especially in the case of high pressures. This means that nonlinear dependence exists between P virial and the volume change ratio ( ∆V V 0 ), which makes no mechanical sense. The nonlinear dependence between P hyd and P virial is caused by the use of constant Ω i . Indeed, the volume of a discrete element should normally decrease with compression. To overcome this problem, a variable volume Ω i is now considered. Ω i is computed as the volume of a sphere having r mean i radius:

r mean i = 1 n v n v ∑ j=1 l ij 2 (5.11)
where l ij = l ij is the length of the beam linking the discrete elements i and j.

Ω i is updated each time these lengths are changed. Figure 5.24 shows that with updated Ω i the dependence between P hyd and P virial becomes practically linear. This dependence can be approximated by a first-order polynomial (5.12).

P hyd = 1 c P virial (5.12)
where c is a parameter, equal to c = 0.725. This parameter can be interpreted as a correction factor of Ω i :

Ω corr i = c Ω i .
As in the case of others microscopic parameters, beyond n p = 10 000, c depends only on the geometric properties (homogeneity and isotropy) of the DEM domain. Hereafter, the corrected-updated expression of Ω i will be used to compute the discrete pressure. In this circumstance, the equivalence between the discrete and continuum pressures is ensured for n p ≥ 10 000.

P i virial = 1 3 trace 1 2 Ω corr i n v ∑ j=1 1 2 l ij ⊗ f ij + f ij ⊗ l ij (5.13)
where Ω corr i is the corrected volume of a sphere having r mean i (5.11) radius. Therefore, no calibration tests are required to apply this model (discrete-continuum pressure equivalence is automatically ensured).

Modeling of nonlinear elasticity

Based on the corrected expression of the virial hydrostatic pressure (5.13) which equals the macroscopic hydrostatic pressure, it is straightforwards to model the nonlinear elastic behavior. For each time step, the virial hydrostatic pressure P i virial at each discrete element i is computed and compared to P nle . If P i virial < P nle , the associated macroscopic properties (corresponding to P = P i virial ) are determined, and hence the microscopic properties (using (5.4) and (5.5)). Finally, the microscopic properties of all the beams connected to i are updated, regardless of their directions.

Modeling of densification

Based on (5.13), a new densification model inspired by the "Newton-Raphson" technique for plasticity is developed (Fig. 5.25). For each time step, the numerical virial hydrostatic pressure (P i virial ) at each discrete element i is computed and compared to the expected one (P ex p ) given experimentally (from the P -∆V/V 0 curve for example). In the case of nonzero difference between these quantities (∆P i = P ex p -P i virial = 0), P i virial must be corrected so that it becomes equal to P ex p . To achieve this aim, it is necessary to derive the appropriate correction forces ∆ f ij that must be injected in the beams linking the discrete elements i and j such that the resultant hydrostatic pressure becomes equal to the expected one. First, ∆P i is divided into three parts:

∆P i = 1 3 (∆σ i xx + ∆σ i yy + ∆σ i zz ) (5.14)
where ∆σ i xx , ∆σ i yy and ∆σ i zz represent the correction amounts of x, y and z directions, respectively. These parts are determined in a way proportionate to their directions contribution for P i virial computation:

∆σ i kk = σ i kk P i virial ∆P i , with k ∈ [x, y, z] (5.15)
The correction amounts of the different space directions are now determined. However, these are given in the form of stresses (useless in the DEM which accepts only forces and torques). The question that arises here is how to translate these correction stresses into correction forces in the beams. In a similar way, the correction stress ∆σ i kk in a direction k is divided into correction stresses of the beams connected to i as follows:

∆σ i kk = n v ∑ j=1 ∆σ ij kk (5.16)
where n v is the number of discrete elements connected to i, ∆σ ij kk correction stress in the direction k of the beam connecting the discrete elements i and j. As with correction stresses of the space directions (∆σ i kk ),∆σ ij kk are also determined in a manner proportionate to their beams (connecting i and j) contribution to compute ∆σ i kk .

∆σ ij kk = σ ij kk n v ∑ j=1 σ ij kk ∆σ i kk = σ ij kk σ i kk ∆σ i kk (5.17)
Knowing (5.17), it is possible to determine ∆ f ij that must be introduced in the beams connecting the discrete elements i and j:

∆ f ij k = 2 Ω corr i ∆σ ij kk l ij k k ∈ [x, y, z] (5.18)
l ij is the oriented length of the beam linking the particles i and j, Ω corr i is the corrected-updated volume associated to i. Keeping in mind Equations (5.15) and (5.17), Equation (5.18) can be simplified as follows:

∆ f ij k = ∆P i P i virial f ij k k ∈ [x, y, z] (5.19)
To check the above calculation, we can compute the corrected hydrostatic pressure taking into account the beam forces corrections and verify if it equals the expected one. The corrected virial stress in a discrete element i is given by:

σ i corr = 1 2 Ω corr i n v ∑ j=1 f ij + ∆ f ij ⊗ l ij + l ij ⊗ f ij + ∆ f ij = 1 2 Ω corr i n v ∑ j=1 f ij ⊗ l ij + l ij ⊗ f ij + 1 2 Ω corr i n v ∑ j=1 ∆ f ij ⊗ l ij + l ij ⊗ ∆ f ij = σ i virial + ∆σ i (5.20)
The corrected hydrostatic pressure at the same discrete element is expressed:

P i corr = P i virial + 1 3 ∑ k∈[x,y,z] 1 2 Ω corr i n v ∑ j=1 ∆ f ij k l ij k (5.21)
Isolate the second term on the right hand side in (5.21):

1 3 ∑ k∈[x,y,z] 1 2 Ω corr i n v ∑ j=1 ∆ f ij k l ij k = 1 6 Ω corr i n v ∑ j=1 ∆ f ij x l ij x + ∆ f ij y l ij y + ∆ f ij z l ij z = 1 6 Ω corr i n v ∑ j=1 2 Ω corr i ∆σ ij xx l ij x l ij x + 2 Ω corr i ∆σ ij yy l ij y l ij y + 2 Ω corr i ∆σ ij zz l ij z l ij z = 1 3 n v ∑ j=1 ∆σ ij xx + ∆σ ij yy + ∆σ ij zz = 1 3 n v ∑ j=1 ∆σ ij xx + n v ∑ j=1 ∆σ ij yy + n v ∑ j=1 ∆σ ij zz = 1 3 ∆σ i xx + ∆σ i yy + ∆σ i zz = ∆P i (5.22)
In conclusion, the corrected hydrostatic pressure is the same that the expected pressure: P i corr = P i virial + ∆P i = P ex p (5.23) Concerning the variation of the Young's modulus E with pressure (Fig. 5.17b), the same technique applied to model the nonlinear elasticity is used. Beyond the densification pressure (P d ), the macroscopic Young's modulus associated to the corrected pressure at each discrete element i is determined. Then, based on (5.5), the associated microscopic Young's modulus E µ can be computed (using (5.5)) and updated for all the beams connected to i (regardless of their directions). Consequently, the densified region keeps its isotropy property. It should be emphasized that only the variation of the Young's modulus with densification is considered in dynamics (no sufficient results are found in the literature concerning the variation of Poisson's ratio with densification).

Application: Plates impact

To validate the new model, the high velocity impact tests of silica glass plates have been reproduced using this model. The correction process is performed with respect to the Hugoniot curve which is approximated by a linear piecewise function (Fig. 5.26). This approximated curve is chosen close to that obtained using the first model with calibrated microscopic parameters, so that the associated results can be compared. To reduce the computation time, the prediction-correction phase is activated only between the densification (P d ) and saturation (P s ) pressures. Figure 5.27 presents a comparison between the numerical results obtained using the two silica glass numerical models (beam-based model and the present predictive-corrective model) for different impact velocities. The results of the predictive-corrective model are very close to those obtained using the beam-based model. Note that using this model the equation of state is practically satisfied at each discrete element. Therefore, the second reason mentioned above to explain the discrepancy between the experimental and numerical values of the maximum particle velocity is not correct. This discrepancy is possibly due to numerical filtering of the HFWs.

Discussion

Using the present predictive-corrective model, only one parameter which is the correction factor c of the local volumes associated to the discrete elements has to be determined by calibration. Once this parameter is determined, the discrete-continuum equivalence is automatically ensured, and then the macroscopic behavior (determined experimentally) can be used as reference to correct the current mechanical state. To compare this model to the first one, an approximated Hugoniot curve close to that obtained using the first model after calibration has been used as reference for the prediction-correction phase. The simulation of the symmetrical impacts of silica glass plates (Fig. (ii) the updating of the microscopic parameters is straightforward; (iii) the updating of the microscopic parameters is performed for all the beams connected to the discrete element in which the macroscopic properties are changed (regardless of their directions), the densified region therefore keeps its isotropy property. Due to these advantages, this model is retained to simulate the LSP processing on silica glass. The poor results of the first model in the region of densification could be enhanced using the experimental Hugoniot curve obtained by Sugiura et al. [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF] as reference for the prediction-correction phase.

Brittle fracture modeling

Naturally, the DEM can deal with multi-fracturing and fragmentation phenomena. Consequently, it is widely used in the literature to simulate tribology [START_REF] Richard | Friction Coefficient as a Macroscopic View of Local Dissipation[END_REF][START_REF] Iordanoff | A discrete element model to investigate sub-surface damage due to surface polishing[END_REF][START_REF] Tan | Discrete element method (DEM) modelling of fracture and damage in the machining process of polycrystalline sic[END_REF], impact [START_REF] Magnier | Numerical simulations of impacts using a discrete element method[END_REF][START_REF] Shiu | Discrete element modelling of missile impacts on a reinforced concrete target[END_REF] or fragmentation [START_REF] D'addetta | From solids to granulates -Discrete element simulations of fracture and fragmentation processes in geomaterials[END_REF][START_REF] Carmona | Fragmentation processes in impact of spheres[END_REF] problems. However, the main challenge encountered in DEM simulations of such applications is the development of a fracture model able to quantitatively fit the fracture mechanisms experimentally observed as reported in the continuum mechanics framework [START_REF] Lawn | Fracture and deformation in brittle solids: A perspective on the issue of scale[END_REF]. As will be seen in this section, even the existent fracture models used in discrete methods give acceptable results at the macroscopic scale, they cannot reproduce the cracking mechanisms at the microscopic scale. To overcome this problem, another fracture model based on the virial stress has been developed in this work. This model is designed for brittle materials (e.g. silica glass) in which completely brittle fracture occurs. Note that this model has been developed in collaboration with D. André [5].

Standard fracture model and its limitations

The approaches most commonly used to model fracture in discrete methods are based on the computation of bond strains [START_REF] Carmona | Fragmentation processes in impact of spheres[END_REF] or stresses [START_REF] Potyondy | Abonded-particle model for rock[END_REF]. Strain and stress can be regarded as dual variables; therefore, these two approaches are fairly similar. The fracture model initially implemented in GranOO (the DEM workbench) is inspired by these ones. It is based on the computation of the Rankine stress which is defined such fracture type is sought to be initiated under tensile stress in mode I. Therefore, the new criterion postulates that a discrete element i is released from its neighbors when the hydrostatic stress computed using (5.13) is positive (tension state in i) and exceeds the fracture strength (σ f ). In this case, all the beam bonds linking this particle to its neighbors break and do not intervene any more to compute the inter-particle forces and torques at the next time step. Here, the fracture strength is that measured at the macroscopic scale σ f . Figure 5.30 shows an illustration of the cracking process using the new criterion.

To validate this criterion, it has been applied to simulate tension, bending and torsion tests using 3D beams. The fracture strength σ f is fixed at 51.5 MPa. To avoid any ambiguity, it should be noted that in our published papers on the fracture model [5,[START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF], this quantity is taken as σ f = 64 MPa. This is because initially we have used constant (not updated) volumes Ω i which leads to σ f = 64 MPa. When using the corrected expression of the virial pressure P virial (5.13), the same quantity becomes σ f ≈ 51.5 MPa. Table 5.5 presents the macroscopic fracture stresses at which there is fracture of the beams for tensile, bending and torsional loading. The numericalmeasured fracture stresses compare favorably with the applied fracture strength. In addition, the crack path for the numerical torsional tests have been observed (Fig.

5.31).

As can be seen, the crack geometry conforms to the prediction by the material strength theory: a crack path is developed along a helical surface oriented at 45 • to the main axis of the cylindrical sample. The above 2D spherical indentation has also been reproduced using this criterion. Figure 5.32 shows the associated result. As expected, the crack pattern exhibits a cone geometry (Fig. 5

.32).

In the previous validation tests of the new fracture criterion, only the elastic behavior has been applied. Now, the complex silica glass behavior is introduced to ensure the proper functioning of this criterion to model fracture of this material; but also to ensure that the silica glass numerical model developed in this work enables to capture the densification effects on the cracking mechanisms. Generally, when a material is indented by a spherical indenter, circular crack surrounding the in- denter appears and propagates to form a cone crack as shown in the 2D spherical indention (Fig. 5.32). However, when a sharp tip (e.g. Vickers indenter) is used, normal glass produces median, radial and lateral cracks. Instead, silica glass has a strong tendency to form a cone crack, even when indented with a sharp tip. Because of its important densification behavior, a spherical densified zone is formed beneath the sharp indenter which in turn operates as a spherical indenter on the rest of the problem domain. A cone crack may therefore sets up which can accompany median, radial and lateral cracks. To study this phenomenon, the Vickers microindentation (Fig. 5.12) has been resumed using high indentation force (30 N). Figure 5.33 presents the numerical cracking response of silica glass indented with a Vickers tip at the beginning of fracture (before the cracking becomes unstable). In this figure, only the discrete elements where the fracture criterion is reached are shown. Qualitatively, the result is in good agreement with those obtained in other studies [7,[START_REF] Gross | Deformation and cracking behavior of glasses indented with diamond tips of various sharpness[END_REF]. Moreover, the radius of the crack obtained numerically (r cone num = 41.14 µm) is in good agreement with the experimental works of Arora et al. [7] (r cone ex p ≈ 40 µm).

Conclusion

To correctly model the silica glass mechanical behavior under high hydrostatic pressure, two models have been proposed for the DEM approach. The first model is based on the normal stress in the cohesive beam bonds. This model involves four microscopic parameters which are the critical nonlinear normal stress σ nle µ , the microscopic densification pressure σ d µ , the microscopic tangential modulus E t µ and the maximum permanent deformation of the beam bonds ε p max µ . These parameters have been determined by calibration. The validation of this model has been performed in statics at both macroscopic and microscopic scales to better analyze the potential problems and avoid the dynamic effects. Relatively good results have been obtained compared with the experimental ones [90,[START_REF] Yoshida | Quantitative evaluation of indentation-induced densification in glass[END_REF]. After, this model has been validated in dynamics by simulation of symmetrical impacts of silica glass plates. In the region of elastic behavior, the numerical results compare favorably with the experimental ones [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF].

Introduction

In the previous chapters, a discrete-continuum coupling approach adapted for dynamic simulations has been developed. This approach couples the Discrete Element Method (DEM) and the Constrained Natural Element Method (CNEM). Subsequently, two numerical models taking into account the various complex phenomena of silica glass mechanical behavior are proposed. These models are adapted for discrete methods. The first model is based on the normal stress in the cohesive beams and involves several microscopic parameters, which have to be determined by calibration tests. The second model is based on the virial stress and is easer to be applied, because only one parameter has to be determined by calibration. Once it is determined, the discrete-continuum equivalence is automatically ensured and the correction of the current mechanical state can be performed with respect to the macroscopic behavior. This makes modeling of complex or even very complex mechanical behaviors relatively simple. Therefore, it is retained in this work to model the silica glass mechanical behavior. On the other hand, a new virial-stress-based fracture model has been proposed in this work to overcome the limitations of the existent fracture models at microscopic scales. This model gives relatively good results at both macroscopic and microscopic scales. All the works detailed so far have been applied to study the laser-induced damage in silica glass through the simulation of Laser Shock Processing on this material. In this application, only a small region is affected by the laser beam and requires a very fine scale analysis. Therefore, the DEM method is only applied in this region; whereas the CNEM is applied beyond it to move away the boundary conditions effects.

A brief description of the LSP test to be simulated is given in Section 6.2. Section 6.3 presents the DEM-CNEM numerical model used to study the silica glass response under laser-induced pressure. The main results of this study are discussed in Section 6.4. Finally, some conclusions and remarks are drawn in the last section.

Brief description of the LSP test to be simulated

A silica glass plate having a square base of 20 × 20 mm 2 and 5 mm thickness is considered in this chapter. To avoid thermal damage of the studied specimen, it is assumed that an aluminum coating layer (opaque layer) is deposited between the silica glass plate and the confining medium (Fig. 6.1). This layer vaporizes during the laser interaction, whose expansion is restricted by the confining medium. The latter (the confining medium) is assumed to be completely transparent with respect to the wavelength of the laser beam used in the process. Then, it is only used to confine the plasma generated by the absorption of the laser beam interacting with the coating layer, so that the peak pressure on the surface of the silica glass can be larger and its duration in time can be longer. Under these circumstances, the thermal damage of the silica glass plate can be neglected and this material can be assumed to undergo only mechanical damage. Also, it is assumed that the laser beam spot diameter is 3 mm and the energy-per-pulse is in the order of several Joules, so that the coating layer is hit by a power density. As mentioned in Chapter 2, several software specialized in the laser-matter interaction can be found. These software are Figure 6.1: Description of LSP test used mainly to evaluate the laser-induced pressure. This part is beyond the scope of this chapter which studies the LSP processing from a mechanical point of view. A temporal pressure profile close to those resulting from these software is used to model the laser-induced pressure.

Numerical model

During the LSP test, only a small region surrounding the laser beam requires a very fine scale analysis. The rest of the plate undergoes only elastic deformation. Therefore, the Discrete Element Method (DEM) is applied only in this region which is of about 6 mm diameter. The rest of the plate is modeled by the Constrained Natural Element Method (CNEM). The two models (DEM and CNEM) overlap in a small region of hollow cylindrical shape with an internal diameter of 4 mm and an external diameter of 6 mm. Figure 6.2 and 6.3 show the associated geometric model. The use of the CNEM around the DEM region intends to move away the boundary condition effects and to prevent spurious reflection of the transversal waves at the lateral border of the DEM region which may distort the longitudinal wave propagation.

The shock surface (surface irradiated by the laser beam) is subjected to a pressure loading, of which the temporal profile is close to those obtained in the literature using laser-matter interaction software. To capture the multiple spalling phenomenon (Fig. 2.14), relatively high pressure with long duration is applied on the shock surface (Fig. 6.4a). The diameter of the shock surface is the same that for the laser beam spot (3 mm). The spatial distribution of pressure on this surface is given in Figure 6.4b, a plateau with gradual decrease at the border to reduce the shear stresses in this zone. To capture the fine scale effects in the shock region, the DEM subdomain is discretized using n p = 50 000 discrete elements (the average radius of the discrete elements is around r mean = 75 µm). If the full studied domain was discretized using Figure 6.2: Geometric coupling model the same r mean , the number of discrete elements would exceed n p = 500 000, which would require several weeks or months of DEM computation just to prepare the DEM domain. So, the coupling approach is advantageous.

To ensure enough CNEM nodes in the overlapping region Ω O , a fine discretization is used for Ω C|Ω O (the part of the CNEM subdomain belonging to Ω O ). The ratio between the discretization characteristic lengths of the DEM and CNEM subdomains in this region is close to 5 (l Ω D c /l Ω C c ≈ 5). Keeping in mind the results of the parametric study (Fig. c . Concerning the junction parameter l, it has been shown that its optimal value l opt is roughly of the same order of magnitude as l Ω O c ; hence, this parameter is fixed at l = l Ω O c . The discrete and continuum energies in the overlapping region are weighted using continuous differential weight functions (Fig. 4.11c). To validate this choice of Arlequin parameters, a preliminary simulation using this geometric model with only elastic behavior has been performed. In this simulation, the high pressure loading is replaced by a moderate pressure of 1GPa to remain within the elastic domain. Figure 6.5 compares the free surface (rear surface) displacements with respect to time obtained with the DEM-CNEM simulation To model the complex mechanical behavior of silica glass in the DEM region, the virial-stress-based model ( §5.4) is used in the DEM region. The correction process in this model is performed with respect to the approximated Hugoniot curve given in Figure 5.26. It should be recalled that using this approximated curve, the virial-stress-based model fails to estimate correctly the plastic wave velocities, compared with those obtained by Sugiura et al. [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF]. This problem can be corrected by using the experimental Hugoniot curve obtained by the same researchers [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF]. However, this curve is given only for moderate pressures (under 20 GPa) . Therefore, the same approximated curve as for the dynamic validation is adopted in this qualitative study. Concerning the CNEM region, only elastic behavior is considered. The fracture of silica glass in the DEM region is modeled by the virial-stress-based fracture model ( §5.5). As seen in Chapter 2, the spalling threshold of silica glass is difficult to measure experimentally and can exceed 2.5 GPa. In this qualitative study, it is fixed at 2.5 GPa.

Results

This section presents the simulation results of the numerical model presented in the previous section. Figure 6.6 shows the pressure waves propagation in the DEM subdomain. The compression pressure wave (blue) is immediately followed by the release wave (red). Arriving at the free surface, the compression pressure wave measures only 8 GPa. Three damping mechanisms can explain this pressure decrease. Two of them are physically-based and introduced in Subsection 2.1.2.4: hydrodynamic damping and elastoplastic damping. The third mechanism is purely numerical in origin and due to filtering of high frequency waves which are not supported by the DEM discretization. Figure 6.7 shows that only an elastic wave has reached the rear surface (no dual wave is observed in the particle velocity profile measured at this surface). The velocity of the elastic wave is about 4961 ms -1 . Moreover, the particle velocity at the rear surface has stabilized at around 117 ms -1 without reversing direction. This means that there is fracture of the material and formation of spalls. The great difference between the maximum and stabilized free surface particle velocities indicates that a thick spall is thought to be formed. Figure 6.8 presents the damage state of the DEM subdomain due to the shock wave propagation. Two consecutive spalls were formed at the rear surface (free surface). In addition, the front surface (shock surface) has also been damaged due to the fast pressure decrease (release wave).

The different results obtained in this section are qualitatively in agreement with the theory of shock wave propagation in materials (presented in Chapter 2) and with the experimental observations on the LSP processing on materials [START_REF] Lescoute | Étude de la fragmentation dynamique de métaux sous choc laser[END_REF][START_REF] Phipps | Laser ablation powered mini-thruster[END_REF].

Conclusion

In this chapter, a qualitative study of the LSP processing on silica glass is given, using the DEM-CNEM coupling approach. To reduce the computation time, only a small region surrounding the laser beam is modeled by the DEM, whereas, the CNEM is applied in the rest of the studied domain to move away the boundary conditions effects. To correctly model the silica glass behavior in the DEM region where complex phenomena can take place, the virial-stress-based model is applied in the region. In addition, the new fracture model (based on the virial stress) is also used in this region to model cracking mechanisms due to the high pressure shock waves propagation. To capture the multiple spalling phenomenon, relatively high pressure with long duration is used in this work. The temporal profile of this loading is close to those obtained using specialized software in laser-matter interaction. The obtained qualitative results are in agreement with the experimental observations on the LSP processing on materials. This is very promising to move on to the quantitative study of this process which would be one of our future works.

Chapter 7

Conclusion and future work

"Science never solves a problem without creating ten more."

George Bernard Shaw

Development and safe use of high-power lasers in various applications require in-depth studies of the origins and sources of laser-induced damage in materials. Despite the current state of the experimental studies, this request is always one of the major concerns encountered in high-power laser applications. Indeed, experimental work is very costly and is not very effective in terms of providing insightful and complete information. This makes numerical support necessary to overcome the experimental difficulties, to reduce the cost of studies and to better characterize the laser-induced damage in materials. The main goal of the present dissertation, Discrete-continuum coupling method for simulation of laser-induced damage in silica glass, was to corroborate the existent studies by developing a robust numerical tool able to correctly predict the mechanical response as well as the damage mechanisms of silica glass subjected to dynamic stresses and particularly to stresses generated by laser irradiation. To achieve the main requests, this dissertation has laid after a general introduction as follows.

Chapter 2 has introduced the main scientific advancements in the different disciplines interacting to accomplish the goals of this work. First, a basic knowledge of the laser-matter interaction theory has been given. The physics behind this theory can be divided into two steps: generation of shock wave by laser ablation and shock wave propagation within the irradiated material. Only the second step has been studied in this work. The temporal pressure profile generated by laser ablation was assumed to be given and it has been directly applied on the shock surface. Second, the complex mechanical behavior of silica glass has been reviewed. In the light of this review, tow numerical models of the silica glass mechanical behavior have been proposed. Finally, a classification of the different numerical methods used to simulate mechanical problems has been given. These methods can be classified into two classes: discrete methods and continuum methods. The first class is naturally more suited to study complex phenomena involving discontinuities, but it is very time-consuming. The second class overcomes this drawback, but is not well adapted to simulate problems requiring a very fine scale analysis. The advantages of these two classes are largely complementary. Therefore, a discrete-continuum coupling approach has been proposed to study the laser-induced damage in silica glass.

Chapter 3 has zoomed in on these two classes to classify the different discrete and continuum numerical methods with respect to our objectives. Based on this classification, the most appropriate discrete and continuum methods have been selected: the Discrete Element Method (DEM) implemented by André et al. [4] and the Constrained Natural Element Method (CNEM) implemented by Yvonnet et al. [START_REF] Yvonnet | Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoi constraint-Méthode des éléments C-Naturels[END_REF].

Chapter 4 has detailed how the selected methods (DEM and CNEM) have been coupled. Based on the Arlequin technique, a 3D DEM-CNEM coupling approach has been developed. Several coupling parameters are involved in this approach. Therefore, to simplify their setting in practice a parametric study of these parameters has been performed. This study has allowed to retain several recommendations simplifying the setting of these parameters in practice. Finally, this new coupling method has been validated using several reference dynamic tests. From an IT point of view, since the GranOO libraries are written in C++ language, the DEM side coupling algorithm has been implemented in the same language. The NESSY libraries are written in C++ language which are interfaced with Python moduli. Therefore, the CNEM side coupling algorithm has been implemented in Python language. The communication between the two codes (GranOO and NESSY) is ensured by a DC-Interface developed in this work in Python language. DC-Interface communicates directly with NESSY (same language) and via an Inter Process Communication (IPC) with GranOO.

Chapter 5 has devoted to the modeling of the silica glass mechanical behavior. To take into account the different specificities of this complex behavior, two models have been proposed. The first model is based on the normal stress in the cohesive beams. The static validation of this model has been performed at both macroscopic and microscopic scales. Relatively good results have been obtained compared with the experimental ones. Also, it is validated in dynamics using symmetrical impacts of silica glass plates. In the region of elasticity, the numerical results compare favorably with the experimental ones. However, relatively poor results have been obtained in the region of densificatiton. This is due to difference between the experimental and calibrated Hugoniot curves. Many microscopic parameters are involved in this model, which makes the calibration process tedious and very time-consuming. To overcome the limitations of this model, a new silica glass numerical model based on the "prediction-correction" technique has been developed. Using this model, only one parameter (the correction factor c of the local volumes associated to the discrete elements) has to be determined by calibration. Once this parameter is determined, the discrete-continuum equivalence is automatically ensured, and then the macroscopic behavior (determined experimentally) can be used as reference to correct the current mechanical state of the simulated problem. To compare this model to the first one, an approximated Hugoniot curve close to that obtained using the first model has been used in this model as reference for the prediction-correction phase. The simulation of the symmetrical impacts of silica glass plates have been reproduced using this new model. The associated results are very close to those obtained using the first model. This new model has been retained to simulate the LSP processing on silica glass. After modeling the silica glass mechanical behavior, a new virial-stressbased fracture model has been proposed in this work to overcome the limitations of the existent fracture models at the microscopic scale. The application of this model to simulate various fracture problems shows its effectiveness at both macroscopic and microscopic scales. All these developments have been implemented in GranOO by adding plug-ins written in C++ language using the advanced Object Oriented (OO) architecture.

Finally, these models have been applied to qualitatively study the laser-induced damage in silica glass through the simulation of Laser Shock Processing on this material using the DEM-CNEM coupling approach. In this application proposed in Chapter 6, the DEM method is applied in a small region surrounding the laser beam; whereas the CNEM is applied beyond this region to move away the boundary conditions effects. This application has demonstrated the interest of the developed coupling method that allows to focus to short scale effects thanks to DEM and to take into account the structural effect thanks to CNEM. The obtained qualitative results are in agreement with the experimental observations on the LSP processing on materials. This is very promising to move on to the quantitative studies of this process.

The different works developed in this dissertation open roads to numerous research directions. A few selected ones are listed bellow:

• Concerning the Laser Shock Processing (LSP), the LSP experiment described in Chapter 6 would be performed. Based on the results of this experiment, the present qualitative studies of the LSP processing on silica glass would be extended to quantitative studies. Also, the extension of the work to study the response of other materials under LSP processing would be helpful in the development and safe use of high-power lasers in different applications.

• The proposed DEM-CNEM coupling approach is a first step toward a longerterm project. The present approach uses the same time step for the two models. As this parameter is controlled by the smallest discretization characteristic length, the CNEM computation would be heavily penalized by using the same time step as for the DEM computation. To overcome this problem, this approach would be extended to use different time step for each model: multitime-step DEM-CNEM coupling approach. Also, in the developed approach, the DEM and CNEM subdomains have to be defined before the computation process. Besides, the DEM model must be applied in all the regions in which fine scale effects can occur. Considering the case of glass scratching simulations [START_REF] Houérou | Surface damage of soda-lime-silica glasses: indentation scratch behavior[END_REF][START_REF] Houérou | Scratchability of Soda-Lime Silica (SLS) glasses: Dynamic fracture analysis[END_REF], the regions requiring a fine scale analysis are relatively large compared to the whole problem domain. Application of the present DEM-CNEM coupling approach would not be advantageous. This has set the groundwork for future research, which should focus on developing a coupling method with moving DEM model to follow the fine scale effects (Fig. 7.1).

• Another direction for future work would be to consider the question of the calibration process of the elastic microscopic parameters in the DEM approach. This process which is very tedious and time-consuming is required each time the modeled material is changed. Therefore, it would be very helpful to develop a mathematical model allowing a simple and straightforward setting of these parameters for a given material.

• The virial-stress-based model of silica glass mechanical behavior gives promising results compared to the beam-based model and would be more investigated at both microscopic and macroscopic scales. Then, the Vickers microindentation of silica glass would be reproduced using this model to enhance the numerical results.

©Figure 1 . 1 :

 11 Figure 1.1: ICF target: a tiny capsule containing atoms of deuterium and tritium is fixed inside a gold cylinder called hohlraum.

  (a) Atmosphere formation: Laser beams or laser-produced Xrays rapidly heat the surface of the fusion target, forming a surrounding plasma envelope. (b) Compression: Fuel is compressed by the rocket-like blowoff of the hot surface material. (c) Ignition: During the final part of the capsule implosion, the fuel core reaches 20 times the density of lead and ignites at 100,000

Figure 1 . 2 :

 12 Figure 1.2: Steps of the inertial confinement fusion (ICF)

Figure 1 . 4 :

 14 Figure 1.4: Global methodology of this dissertation

Figure 2 . 1 :

 21 Figure 2.1: Electromagnetic spectrum

Figure 2 . 2 :Figure 2 . 3 :

 2223 Figure 2.2: Physical processes of laser matter interaction

Figure 2 . 4 :

 24 Figure 2.4: Physical processes of material laser ablation

Figure 2 . 5 :

 25 Figure 2.5: Evolution of a pressure pulse during its propagation in a material

Figure 2 . 7 :

 27 Figure 2.7: Experimentally measured EOS curve of V s vs. V p for several materials. (Taken from [126])

Figure 2 . 9 :

 29 Figure 2.9: Permanent deformation effects on the Hugoniot and the wave propagation; left: Hugoniot curve (solid line), isentropic hydrostatic compression (dashed line); right: shock wave instability due to the HEL presence (dual wave).

Figure 2 . 10 :

 210 Figure 2.10: Illustration of the phenomena spoiling the propagation of shock waves

  (a) Phenomenological description of the spalling due to the propagation of a shock compression (b) Multiple spalling process

Figure 2 .Figure 2 . 14 :

 2214 Figure 2.13: Laser-induced spalling

Figure 2 . 15 :

 215 Figure 2.15: Relationships between the elastic moduli of silica glass and the hydrostatic pressure; K, G and C 11 are respectively the bulk, shear and longitudinal moduli. (Taken from [95])

  2.16).

Figure 2 . 16 :

 216 Figure 2.16: Variations in mechanical properties of silica glass with densification (Experimental works of Ji [90])

Figure 2 . 18 :Figure 2 . 19 :

 218219 Figure 2.18: Particle velocity profiles in silica glass under shock compression: nonlinear elastic behavior, linear elastic behavior and apparition of a second wave beyond the HEL. (Taken from [168])

Figure 2 . 20 :Figure 2 . 21 :

 220221 Figure 2.20: Hugoniot curve of silica glass obtained experimentally (Taken from [129])

Figure 2 . 22 :

 222 Figure 2.22: Different steps of numerical simulation

Figure 2 . 23 :

 223 Figure 2.23: Example of cohesive links used in DMs

Figure 2 . 24 :

 224 Figure 2.24: DEM modeling of a continuum

2 . 23 )

 223 can be of different natures according to the physical properties of the material being modeled. Figure2.24 illustrates an example of a continuum domain modeled with DEM. Application of DMs in continuum media simulations faces two significant challenges:

Figure 2 . 25 :

 225 Figure 2.25: Types of coupling methods

Figure 3 . 2 :Figure 3 . 3 :

 3233 Figure 3.2: Classification of continuum methods

Figure 3 . 4 :

 34 Figure 3.4: Eulerian mesh for a rabbit simulation

Figure 3 . 5 :

 35 Figure 3.5: Meshfree discretization for rabbit simulation

Figure 3 . 6 :

 36 Figure 3.6: Different geometries of problem domains

Figure 4 . 1 :

 41 Figure 4.1: Illustration of DEM modeling of a continuum

4 . 2 :

 42 Free random filling of the DEM domain during the first step

Figure 4 . 3 :Figure 4 . 4 :

 4344 Figure 4.3: Geometrical arrangements of a 3D sphere packing with different values of X (Taken from [3])

Figure 4 . 5 :

 45 Figure 4.5: 3D histograms for different DEM domains obtained using different radial dispersions X (Taken from [3])

Figure 4 . 6 :

 46 Figure 4.6: Cohesive beam bond between two discrete elements (Taken from [4])

Figure 4 .

 4 7 presents an example of Voronoï diagram and the associated Delaunay tessellation in 2D.

Figure 4 . 7 :

 47 Figure 4.7: Geometric construction of the Voronoï diagram and the associated Delaunay tesselation for a set of nodes in the plane (2D)

  where n v is the number of the natural neighbors of x. The left view of Figure4.8 gives and illustration of 2D calculation of Sibson shape functions. In 2D, the Sibson shape function φ i associtated to a node i is calculated in x as the ratio of the A and B polygon areas (Fig.4.8):

Figure 4 .

 4 [START_REF] Aubertin | A coupled molecular dynamics and extended finite element method for dynamic crack propagation[END_REF] presents an example of a 2D NN shape function support. Contrary to other meshfree methods, the support domains of the NN shape functions are automatically defined by the natural neighbors. These supports are always

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Calculation of Sibson shape functions (left) and Laplace shape functions (right)

Figure 4 . 10 :

 410 Figure 4.10: Global problem domain decomposition

Figure 4 . 11 :

 411 Figure 4.11: Examples of weight functions

Figure 4 .

 4 [START_REF] Bauman | On the application of the Arlequin Method to the Coupling of Particle and Continuum Models[END_REF] presents several examples of these functions.

Figure 4 . 12 :

 412 Figure 4.12: Discrete-Continuum decomposition

(

Figure 4 . 13 :

 413 Figure 4.13: Different configurations of the discretized subdomains; blue squares: CNEM nodes; red circles: DEM particles; black crosses: nodes of the overlapping subdomain

  {λ} is the vector of coupling forces. C L2 C and C H1 C are the continuum L2 and H1 coupling matrices.

54 )

 54 Solving Equation (4.52), {λ} n+1 can be obtained.• Computation of F λ n+1 and f λ n+1 using Equation (4.50).• Computation of the linear velocities .

Figure 4 .

 4 [START_REF] Belikov | The non-Sibsonian interpolation : a new method of interpolation of the values of a function on an arbitrary set of points[END_REF] gives the general architecture of the coupling approach.

Figure 4 . 14 :

 414 Figure 4.14: Coupling architecture

4 . 16 ).

 416 As shown in Figure 4.16b, the Fourier spectrum contains powerful high frequency waves (greater than f CNEM c

Figure 4 . 15 :

 415 Figure 4.15: Beam model for the parametric study

Figure 4 . 16 :

 416 Figure 4.16: Tensile loading used for the parametric study and the associated spectral analysis (computed from FFT)

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

 418434521 Figure 4.18: Conditioning of [A] with respect to l; L O = 6 mm, coarse mediator space, continuous weight functions with ε = 0.05

Figure 4 . 19 : 5 Figure 4 . 20 : 5 4. 5 . 2 . 2

 41954205522 Figure 4.19: Free-end displacements obtained using DEM, CNEM and the coupling method, and the associated spectral analyses (computed from FFT); L O = 2 mm, fine mediator space, constant weight functions α = 0.5

Figure 4 . 21 :

 421 Figure 4.21: Free-end displacements obtained using the coupling method with different constant weight functions: α = 0.3, α = 0.5 and α = 0.8, and the associated spectral analyses (computed from FFT); L O = 2 mm, fine mediator space

Figure 4 . 22 :

 422 Figure 4.22: Linear velocities at the check points for different constant weight functions: α = 0.3, α = 0.5 and α = 0.8 ; L O = 2 mm, fine mediator space

Figure 4 . 23 :Table 4 . 4 :

 42344 Figure 4.23: Free-end displacements obtained using the coupling method, and the associated spectral analyses (computed from FFT) for different values of ε; L O = 2 mm, fine mediator space, continuous weight functions

4 . 11 )

 411 to reduce the weight jump in the vicinity of ∂Ω O (boundary of the overlapping zone). Figures 4.33

  , 4.34 and 4.35

Figure 4 . 24 :

 424 Figure 4.24: Free-end displacements obtained using DEM, CNEM and the coupling method, and the associated spectral analyses (computed from FFT); L O = 2 mm, fine mediator space, continuous weight functions

Figure 4 .

 4 Figure 4.25: (a) Linear velocities at the different check points for the case of continuous α; (b) comparison of linear velocities at the "DemCheckPoint" obtained using continuous and constant (α = 0.5) weight functions; L O = 2 mm, fine mediator space

Figure 4 . 26 :

 426 Figure 4.26: Velocities comparison in the overlapping zone using fine and coarse mediator (multiplier) spaces ; L O = 2 mm, continuous weight function, ε = 0.05

  and 4.28 provides the coupling results for different L O using fine and coarse mediator spaces, respectively. In the case of fine M h

4. 5 . 5

 55 Dependence between L O and M h O This subsection highlights the dependence between L O and M h O . First, the discretization characteristic length of the overlapping subdomain l Ω O c will be investigated. The width of the overlapping region is fixed at L O = 2 mm, the discretization characteristic lengths l Ω C c and l Ω D c of respectively the continuum and discrete subdomains are chosen as: l Ω C c = 5 l Ω D c . Figure 4.29 shows the free-end displacements for different choice of l Ω O c . According to this figure, the optimal value of l

Figure 4 . 27 :Figure 4 . 28 :

 427428 Figure 4.27: Free-end displacements obtained using the coupling method for L O = 2 mm, L O = 4 mm and L O = 6 mm.; Fine mediator space , continuous weight functions, ε = 0.05

Figure 4 .

 4 31 presents the associated results.

Figure 4 . 29 :Figure 4 . 30 :Figure 4 . 31 :

 429430431 Figure 4.29: Influence of the discretization characteristic length of the overlapping subdomain l Ω O c , L O = 2 mm, l Ω C c = 5 l Ω D c

Figure 4 .

 4 Figure 4.32: Validation model

  4.33, 4.34, 4.36 and 4.35 and Tab. 4.5) validates the new coupling method.

Figure 4 . 33 :

 433 Figure 4.33: Tensile loading: Free-end displacements U x mean obtained by DEM, CNEM and the coupling method; L O = 10 mm, coarse mediator space , continuous differentiable weighting, ε = 0.05

Figure 4 . 34 :

 434 Figure 4.34: Bending loading: Free-end displacements U y mean obtained by DEM, CNEM and the coupling method; L O = 10 mm, coarse mediator space , continuous differentiable weighting, ε = 0.05

Figure 4 . 35 :Figure 4 . 36 :

 435436 Figure 4.35: Torsional loading: Free-end rotation θ x mean obtained by DEM, CNEM and the coupling method; L O = 10 mm, coarse mediator space , continuous differentiable weight functions, ε = 0.05

Figure 5 . 2 :Figure 5 . 3 :

 5253 Figure 5.2: Rheological model of the cohesive bonds used to model densification of silica glass

Figure 5 . 4 :Figure 5 . 5 :

 5455 Figure 5.4: Linear piecewise (C 0 ) approximation of the variations of the silica glass mechanical properties with compressive pressure

Figure 5 . 7 :

 57 Figure 5.7: Sensitivity of the macroscopic densification behavior to the number of discrete elements (n p )

σ d µ = 5 Figure 5 . 8 :σ d µ = 5 Figure 5 . 9 :Figure 5 . 11 :

 558559511 Figure 5.8: Variation of macroscopic densification pressure (P d ) with the microscopic densification parameters: microscopic densification pressure (σ d µ ), microscopic tangential modulus (E t µ ) and maximum permanent deformation of the beam bonds (ε p max µ )

5 . 11 ).

 511 As can be observed in Figure5.11, the numerical results are in agreement with the experimental ones.Microscopic validation: Vickers micro-indentationSilica glass experiences permanent densification under high compressive pressure. Comparable stresses can be reached beneath a sharp indenter (e.g. Vickers indenter) during indentation test. To verify the validity of the developed densification model at microscopic scale, Vickers micro-indentations are simulated. Low indentation forces are used to study the densification behavior in the absence of fracture: F indent = 0.1 N and F indent = 0.5 N.

Figure 5 . 14 :

 514 Figure 5.14: Top view of indentation print (blue: negative volume change) and piledup material around indenter (red: positive volume change). Left: indentation force F indent = 0.1 N; right: indentation force F indent = 0.5 N

  Figure 5.15: Cut-away view of the indented region; F indent = 0.5 N; (a) experimental results from [90], (b) DEM numerical results; (c) DEM-CNEM coupling numerical result

19

 19

  and 2.20, which show the variation of the bulk modulus with the density change ( ∆ρ ρ 0 ) and the variation of the longitudinal stress with density change (Hugoniot curve).

Figure 5 . 18 :Figure 5 . 19 :

 518519 Figure 5.18: Hugoniot curve: silica glass numerical model based on the normal stress in the cohesive beams

( a )Figure 5 . 20 :

 a520 Figure 5.20: Impact velocity loading

Figure 5 . 22 :roposal 1 )Figure 5 . 23 :

 5221523 Figure 5.22: Particle velocity profiles corresponding to impact velocities of 880 ms -1(black), 2580 ms -1 (red) and 3310 ms -1 (green): comparison between numerical and experimental results.

Figure 5 . 24 :

 524 Figure 5.24: Variation of virial hydrostatic pressure with macroscopic hydrostatic pressure for constant and updated discrete element volume (Ω i )

Figure 5 . 25 :

 525 Figure 5.25: Predictive-Corrective model

Figure 5 . 26 :

 526 Figure 5.26: Hugoniot curve approximation for the predictive-corrective densification model

( a )

 a Impact velocity: V Impact = 880 ms -1 (b) Impact velocity: V Impact = 2580 ms -1 (c) Impact velocity: V Impact = 3310 ms -1

Figure 5 . 27 :

 527 Figure 5.27: Particle velocity profiles in silica glass: comparison between the numerical results of the beam-based model and the predictive-corrective model.

Figure 5 . 30 :Table 5 . 5 :

 53055 Figure 5.30: The cracking mechanisms at the microscopic scale using the new fracture model

Figure 5 . 31 :

 531 Figure 5.31: View of crack path in a torsional test; the discrete elements in which the fracture criterion is fulfilled are highlighted.

Figure 5 . 32 :

 532 Figure 5.32: Qualitative 2D indentation test showing cracks produced when using the new fracture criterion.

Figure 5 . 33 :

 533 Figure 5.33: Cracking pattern of silica glass indented with a Vickers tip; F indent = 30 N; only the discrete elements where the fracture criterion is reached are shown in (b) and (c)

4 . 29 )

 429 , the discretization characteristic length of the overlapping subdomain lΩ O c is chosen as l Ω O c = 2 l Ω D

Figure 6 . 3 :Figure 6 . 4 :

 6364 Figure 6.3: 3D visualization of the different coupling subdomains

Figure 6 . 5 :

 65 Figure 6.5: Time evolution of the free surface (rear surface) displacements obtained using the DEM-CNEM and CNEM approaches.

8 Figure 6 . 6 :Figure 6 . 8 : 1 )Figure 6 . 7 :

 86668167 Figure 6.6: Pressure waves propagation in the DEM region; blue: compression wave; red: release wave

Figure 7 . 1 :

 71 Figure 7.1: Simulation of silica glass scratching using adaptive DEM model

  2.1) extends from low frequencies used for modern radio communication to high frequencies encountered in nuclear decay, e.g. gamma radiation. The high-power lasers which are the concern of this dissertation emit radiations between infrared and ultraviolet.

	Gamma ray	X-ray	Ultraviolet	Visible	Infrared	Microwave	Radio

.19). Figures 2.20 and 2.21 present respectively the Hugoniot curve and shock polar of

  silica glass under shock compression.

  O i are the CNN shape functions defined on D ad,h , U ad,h and M h , respectively. d i are the DEM particle displacements, u i are the CNEM nodal displacements and λ i are the nodal Lagrange multipliers.

	37)
	where n D , n C and n O are respectively the total number of nodes in Ω h D|Ω O , Ω h C and Ω h O . φ D i , φ C i and φ

Using (4.35)

,

(4.36)

and (4.37), the discretized forms of the different weak formulations (4.32), (4.33) and (4.34) can be derived. These discretized equations can

  CNEM code {λ} n+1 ← Resolution of system (4.52) {F c } n+1 ← Application of Eq. (4.50) { f c } n+1 ← Application of Eq. (4.50) { u} n+1 ← Application of Eq. (4.49) ḋ n+1 ← Application of Eq. (4.49)

	55)
	• Transfer of { ü} n+1 and . d n+1 to the DEM code. . u n+1 to the CNEM code and transfer of d n+1 and
	Algorithm 3 Interface Algorithm
	Require: Coupling matrices
	n ← 0 for all iteration n do ḋ * n+1 , d * n+1 ← DEM code { u} *
	{ ü} n+1 ← Application of Eq. (4.55) d n+1 ← Application of Eq. (4.55) ḋ n+1 , d

n+1 , { ü} * n+1 ← n+1 , → DEM code

Table 4 . 2

 42

	Hz)

: Comparison of DEM, CNEM and analytical results is better than the L2 coupling (4.29). However, contrary to what is presented in the literature, the H1 coupling becomes worse if l exceeds some small value.

Table 4 .

 4 3 presents the optimal values of the junction parameter (l opt ) determined numerically for different L O and l Ω O c

Table 5 .

 5 1 gives the calibrated microscopic densification parameters that must be used to ensure the correct macroscopic densification behavior of silica glass.

	σ d µ (GPa) E t µ (GPa) ε 7.13 98.82 0.067 p max µ
	Table 5.1: Static microscopic parameters of silica glass mechanical behavior (obtained
	by calibration)

Table 5 .

 5 2: Vickers indentation results: low indentation forces

Table 5 .

 5 4: Shock fronts velocities for various impact velocities

			st shock front	2 nd shock front
			Shock	Particle	Shock	Particle
			velocity	velocity	velocity	velocity
			( ms -1 )	( ms -1 )	( ms -1 )	( ms -1 )
	880	Num Exp [168]	4931 4940	434 440	--	--
	2580	Num Exp [168]	4884 5030	802 750	4230 3910	1275 1290
	3310	Num Exp [168]	4901 4990	806 800	4300 4550	1630 1655

  5.19) have been reproduced using this new model. The associated results are very close to those obtained using the first model. Compared to the beam-based model (first model), the present model has several advantages: (i) since only one parameter have to determined by calibration, this model is easy to be applied;

http://www.granoo.org

http://plateformesn-m2p.ensam.eu/SphinxDoc/cnem/cnem3d.html

(a) Small number of discrete elements n p (b) Great number of discrete elements n p

Acknowledgements

First and foremost, I would like to thank my supervisors Ivan IORDANOFF, Jean-Luc CHARLES and Frédéric DAU. They have been a constant source of support all along the last three years.

• To start with Ivan, I admire the way he approaches a problem, his initial preprocessing by removing all the "garbage". This has been of crucial help in the problem solving involved in the thesis.

Chapter 3

Choice of the numerical methods

"The choices we make by accident are just as important as the choices we make by design."

Shad Helmstetter

Abstract

This chapter aims to select the appropriate discrete and continuum methods to be coupled in order to simulate the laser-induced damage in silica glass. First, an overview of the most commonly used numerical methods in Computational Mechanics will be given. Then, these methods will be classified according to their advantages and drawbacks with respect to our goal. Using this classification, the methods that best meet the expectations of the present work will be chosen. "" Figure 5.1: Variations of the silica glass mechanical properties in the region of nonlinear elasticity (Deduced from experimental works of Kondo et al. [START_REF] Kondo | Nonlinear pressure dependence of the elastic moduli of fused quartz up to 3 GPa[END_REF]) function of rµ for a constant E µ .

According to these facts, the relationships between the microscopic and macroscopic parameters can be written as in (5.4) and (5.5): rµ = a 0 + a 1 ν + a 2 ν 2 (5.4) The various parameters of Equations (5.4) and (5.5) are numerically determined using the Moving Least Squares (MLS) technique. Bearing in mind the simplifying assumption (5.2), the variation of the macroscopic properties with pressure translates at the level of cohesive beams into variation of the microscopic properties with normal stress in the beams. Based on Equations (5.4) and (5.5), the variation of the microscopic properties E µ and rµ in the region of nonlinear elasticity can be determined. At the level of cohesive beams, the region of nonlinear elasticity is expressed in terms of normal stress as:

, where σ nle µ is the normal stress threshold from which the elastic behavior becomes linear. In practical terms, the variations of the macroscopic properties with pressure in the region of nonlinear elasticity are expressed in terms of dimensionless pressure P P nle . The microscopic parameters of a beam at a given σ n µ (σ n µ < σ nle µ ) correspond to the macroscopic properties at P, such that

For each beam, the normal stress σ n µ is computed and compared to σ nle µ each time step. If σ n µ < σ nle µ , the microscopic properties of this beam must be updated: first, the associated macroscopic properties which correspond to P = σ µ σ nle µ P nle must be determined; then, the new microscopic properties can be evaluated, using (5.4) and (5.5), and updated. Here, the critical microscopic value of the The sphere model was discretized using various numbers of discrete elements n p (from 100 to 20 000). The microscopic densification parameters were fixed as follows: the microscopic densification pressure σ d µ = 5 GPa, the microscopic tangential modulus E t µ = 0 GPa and the maximum permanent deformation of the beam bonds ε p max µ = 0.05. Figure 5.7 presents the results of this study. For a small number of discrete elements (n p ), the densification response greatly fluctuates with this parameter (Fig. 5.7a). These fluctuations can be reduced by increasing n p . Beyond n p = 10 000, the densification response becomes very weakly affected by the discrete element number (Fig. 5.7b). This value of n p is in agreement with the results of André et al. [3,4] who have shown that beyond n p = 10 000, the homogeneity and isotropy of the discrete domain (DEM domain) are ensured and the elastic behavior becomes nearly independent of n p .

In the remainder of this study, n p = 10 000 discrete elements will be used to discretize the numerical samples. Figure 5.8a shows that the macroscopic densification pressure P d depends linearly on the microscopic one σ d µ . In addition, as observed in Figures 5.8b and 5.8c, P d is effectively independent of the other microscopic properties, i.e. E t µ and ε p max µ . Figure 5.9 shows the influence of the various microscopic parameters on the macroscopic saturation pressure P s . This parameter (P s ) depends linearly on the microscopic densification pressure σ d µ , also it slightly depends on the microscopic tangential modulus E t µ (around 12% P s variation for E t µ up to 100 GPa). However, it is effectively independent of the maximum permanent deformation of the beam bonds ε p max µ (Fig. 5.9c), only 5% P s variation for ε p max µ up to 0.12. Concerning the third macroscopic parameter ∆V V 0 , Figure 5.10 shows that this parameter depends only on E t µ and ε p max µ , and it is effectively independent of σ d µ . After studying the influence of the microscopic parameters on the macroscopic densification behavior, calibration becomes fairly easy. Because the macroscopic densification pressure P d depends only on the microscopic one σ d µ , this parameter must be calibrated first. Subsequently, the microscopic tangent modulus E t µ can be calibrated to ensure the expected value of the macroscopic saturation pressure. Finally, the maximum permanent deformation of the beam bonds ε p max µ can be set to find the appropriate global as:

where σ R µ is the maximum equivalent Rankine stress, σ n max µ is the maximum normal stress and and τ max µ is the maximum shear stress. If the maximum equivalent Rankine stress of a cohesive beam σ R µ is larger than a calibrated critical value σ f µ , the considered beam is deleted. This model has been investigated using a qualitative 2D indentation test with a spherical indenter [5]. Figure 5.28 shows the associated qualitative results using the standard fracture criterion defined by (5.24). This simulation does not produce the hertzian cone crack as experimentally observed (Fig. 5

.29).

Several cracks initiate beneath the indenter and propagate vertically throughout the thickness of the material being simulated. In conclusion, even this criterion is moderately satisfactory [3,[START_REF] Srivastava | Strength of Materials[END_REF] at the macroscopic scale, it cannot reproduce correctly the cracking mechanisms at the microscopic scale.

Fracture model based on the virial stress

Based on the virial stress, a new fracture model has been developed. Remember that in this work it is assumed that completely brittle fracture occurs in silica glass, and However, relatively poor results have been obtained in the region of densification. In addition, this model presents some drawbacks. Many microscopic parameters are involved in this model, which makes the calibration process tedious and very timeconsuming. Using this model, only the beams in which σ d µ is reached experience the variation of the microscopic properties with normal stress. A state of mechanical anisotropy can therefore be established in the material being simulated.

To overcome the limitations of this model, another silica glass numerical model has been developed. This model is based on "prediction-correction" technique, using the virial stress at the discrete elements. Only one parameter, which is the correction factor c of the local volumes associated to the discrete elements, has to be determined by calibration. Once this parameter is determined, the discrete-continuum equivalence is automatically ensured. Therefore, the macroscopic behavior (determined experimentally) can be used as reference to correct the current mechanical state. To compare this model to the first one, an approximated Hugoniot curve close to that obtained using the first one after calibration is used as reference for the predictioncorrection phase. The simulations of the symmetrical impacts of silica glass plates have been reproduced using this new model. The associated results are very close to those obtained using the first model. The poor results of the first model in the region of densification could be improved using the experimental Hugoniot curve obtained by Sugiura et al. [START_REF] Sugiura | Dynamic response of fused quartz in the permanent densification region[END_REF] as reference for the prediction-correction phase.

Compared with the beam-based model, the virial-stress-based model has several advantages: (i) since only one parameter should be calibrated, this model is easy to be applied; (ii) the updating of the microscopic parameters is preformed for all the beams connected to the discrete element in which the macroscopic parameters are changed, the densified region keeps therefore its isotropy property. Due to these advantages, this model is retained to simulate the LSP processing on silica glass. Also, to circumvent the limitations of the existent fracture models at the microscopic scale, a new virial-stress-based fracture model has been developed in this chapter. The application of this model on various fracture problems shows its effectiveness at both macroscopic and microscopic scales.

All the developments seen in this chapter have been implemented in the GranOO workbench by adding plug-ins written in C++ language using the advanced Object Oriented (OO) architecture.

Chapter 6

Simulation of Laser Shock Processing on silica glass: qualitative study "No amount of experimentation can ever prove me right; a single experiment can prove me wrong."

Albert Einstein

Abstract

The present chapter aims to study the laser-induced damage in silica glass through simulation of the Laser Shock Processing on this material, using the discrete-continuum coupling approach (DEM-CNEM). Only the mechanical part of this process is studied here. The laser-matter interaction part is beyond the scope of this work. Due to lack of appropriate experimental results, only a quantitative study is proposed in this chapter.