Michel Bellet

Charles-André Gandin

Keywords: Additive manufacturing, modeling, level set, finite elements, fluid dynamics, solid mechanics

This three-year PhD is an unforgettable experience in my life. Cemef provided me a good research environment and opportunities to learn from other colleagues. It is also my honor to work in such a laboratory with good reputation and famous software products in the domain of material forming, both in academic committee and industry.

First of all, I

VP Viscoplastic

κ t Total curvature m -1 λ Thermal conductivity W • m -1 • K -1 λ i , µ i , i = 1,
x ij = x j -x i Edge vector m A Absorptivity c Light speed m • s -1 C p Specific heat capacity J • kg -1 • K -1
q imp Imposed heat flux W • m -2 R Reflection R g Gas constant J • mol -1 • K -1

Context

Recent developments of Additive Manufacturing (AM) technologies opens the door of new manufacturing modes. Among them, the advantages of Selective Laser Melting (SLM) are well recognized. However, most published researches focus on the application to metals while less attention is paid to its application to ceramics. Regardless of the good mechanical performance of ceramics, it is difficult to process them due to their low absorption to commonly used laser like Yb:YAG and their relatively weak resistance to thermal shocks. However, recent studies have clearly identified the possibility to produce dense and crack-free ceramic parts by the SLM technique [START_REF] Hagedorn | Net shaped high performance oxide ceramic parts by selective laser melting[END_REF]. As a result, the project C ÉFAL É (Fabrication additive de pièces C Éramiques hautes performances par Fusion lAser séLEctive) has been proposed with two principal objectives:

• production of part with complex geometry and good mechanical properties without post-treatments;

• design of fine microstructures in order to enforce mechanical properties of ceramics.

This project was originally oriented to the materials of the MAX * phases and the eutectic of Al 2 O 3 -ZrO 2 . However, due to the complexity of these two types of materials, Al 2 O 3 was selected for the essential of the study. The project has been developed as a collaboration between three laboratories. The powder was synthesized by Centre des Sciences des Matériaux et des Structures (Ecole Nationale Supérieure des MINES de Saint-Etienne). Centre des Matériaux (CdM, MINES ParisTech) took the charge of the experimental development of the SLM process, including the characterization of powder and the manufacturing of simple parts with desired microstructures and good mechanical properties. Numerical modeling and simulation were carried out at Centre de Mise en Forme des Matériaux (Cemef, MINES ParisTech) based on welding simulation models previously developed by Hamide [2], Desmaison [START_REF] Desmaison | Modélisation Numérique d'un Procédé de Soudage Hybride Arc / Laser En Approche Level Set : Application Au Soudage Multi-Passes de Tôles d'acier de Forte Épaisseur[END_REF] and Chen [START_REF] Chen | Three Dimensional Cellular Automaton-Finite Element (CAFE) Modeling for the Grain Structures Development in Gas Tungsten / Metal Arc Welding Processes[END_REF] within the numerical library Cimlib at Cemef. The comparison between experiments and numerical simulations was realized on the bases of data issued from experimental works conducted at CdM. The objectives of the numerical modeling and simulation are:

• modeling of SLM process at the scale of track development, including laser/material interaction, heat transfer, fluid dynamics in the melt pool and solid mechanics;

• prediction of temperature field, liquid velocity in the melt pool and stress distribution around the melted zone;

Résumé

Dans ce chapitre, l'histoire et le principe de la fabrication additive sont brièvement présentés , ainsi que leurs avantages et inconvénients rencontrés dans la production. Différents procédés de fabrication additive sont comparés en fonction de leurs différents principes et applications. Le procédé SLM est ensuite détaillé par les étapes du procédé, le laser appliqué, les matériaux et les applications.

Les multiples phénomènes physiques apparaissant dans SLM sont présentés. Les phénomènes de transfert de chaleur par conduction, convection, rayonnement et évaporation sont détaillés. Le transfert de masse provoqué par l'éjection du matériau et les éclaboussures de liquide est considéré comme ayant un effet important sur la morphologie du cordon. Concernant la dynamique du bain liquide, la tension de surface, l'effet de Marangoni et la pression de recul exercent une influence significative sur la forme du bain liquide et aussi la morphologie du cordon. La formation de fissures est liée au gradient thermique élevé dans le procédé SLM et la méthode d'élimination par préchauffage est présentée dans la littérature. La structure fine due à la solidification rapide et la texture liée à la solidification directionnelle sont également des caractérisations spécifiques du SLM. A la fin, les méthodes de modélisation du procédé SLM aux différentes échelles sont comparées et les challenges de modélisation sont soulignés.

ADDITIVE MANUFACTURING

Additive Manufacturing 1.Principle, advantages and challenges

According to Wohlers and Caffrey, Additive Manufacturing (AM, or 3D printing) is defined as a process of joining materials from 3D model data, usually layer upon layer, as opposed to subtractive manufacturing methodologies [5], such as machining technologies. It began with a rapid prototyping process called stereolithography [START_REF] Hull | Apparatus for Production of Three-Dimensional Objects by Stereolithography[END_REF] in 1980s in the United States. This process uses ultraviolet sensitive liquid polymer to construct a part by layer solidification. From that time, many researches have been carried out and new processes have been developed. However, all processes are based on the same principle -layer by layer construction.

A typical AM process cycle begins with a 3D Computer Aided Design (CAD) model, which can be directly modeled in CAD software or by reverse engineering equipment (e.g. intraoral scanner in dental prostheses [START_REF] Gebhardt | Additive Manufacturing by selective laser melting the realizer desktop machine and its application for the dental industry[END_REF]). The 3D model can be sliced into a bunch of sheets. Each sheet corresponds to a construction layer in AM. The scanning trajectory is then generated to pilot the AM machine and construct the final part. At last, post-processing like surface and heat treatment is often conducted in order to reduce the roughness and optimize the mechanical performance. Compared with conventional manufacturing processes, AM has several advantages [START_REF] Attaran | The rise of 3-D printing: The advantages of additive manufacturing over traditional manufacturing[END_REF]:

• Less limits for the complexity of geometry. Internal channels, lattices or assembled parts can be directly achieved;

• Product variety is easy as one just needs to change the 3D CAD model;

• No time gap between design and prototyping;

• Less waste as material is added rather than subtracted from a rough near-net shape part and the non melted powder may be recycled for powder-bed based processes.

However, there are also some challenges in additive manufacturing [START_REF] Attaran | The rise of 3-D printing: The advantages of additive manufacturing over traditional manufacturing[END_REF]. The first difficulty raised is the low production rate. The construction of part is very time consuming compared with conventional processes like casting and machining. Consequently, additive manufacturing is more adapted to unique parts or to small series. On the other hand, it needs a lot of efforts to process individual parts, including the design of trajectory and construction strategy, the adjustment of process parameters and so on. Also, the poor control of process often ends up with defects like cracking, poor surface quality and low mechanical properties. All these effects lead to high price of final parts. Another restriction is the construction size. For example, most of the commercialized machines have a build chamber dimension limited to about 300 mm. Even though some companies provide machines with larger construction size, the distortion of such large parts remains a big problem. Despite of these disadvantages, additive manufacturing is still an attractive technology, especially due to its capacity to produce parts that may not be made by conventional processes, such as parts with complex geometry like internal cooling channels and lattice structures.

Different processes

The main differences between AM processes are layer deposition technology and the choice of materials. Some principal processes are detailed in the following and corresponding application cases are shown in Figure 1.1.

• Stereolithography Apparatus (SLA)

The research of this process began at the 1970s. It was patented in 1986 [START_REF] Hull | Apparatus for Production of Three-Dimensional Objects by Stereolithography[END_REF] and commercialized by the first 3D printing company, 3D System Inc. The applied material is dedicated to photopolymers. By focusing an ultraviolet laser to a vat of photopolymer resin, the latter is solidified and forms a single layer of desired 3D object [START_REF] Fouassier | Photopolymerization Reactions[END_REF]. A new process based on this technology is the Continuous Liquid Interface Production (CLIP), invented by J. DeSimone et al. [START_REF] Tumbleston | Continuous liquid interface production of 3D objects[END_REF]. According to the authors, CLIP is 100 times faster than traditional SLA. • Fused Deposition Modeling (FDM) This process was developed by S. Scott Crump in the late 1980s and commercialized by Stratasys in 1991 [START_REF] Wohlers | History of additive manufacturing[END_REF]. The part is fabricated by extruding molten material from a nozzle to deposit a layer that solidifies on top of the previous layer. Variable materials can be used in this process, such as acrylonitrile butadiene styrene, polylactic acid, polycarbonate, polyamide, polystyrene, etc. This process is clean and simple to use, usually applied to rapid prototyping. In addition, there is a large open-source development community based on this technology (called RepRap). The low price of machines 1.2. SELECTIVE LASER MELTING makes it accessible to personal project.

• Direct Metal Deposition (DMD) DMD was firstly announced by Precision Optical Manufacturing in 2000 [START_REF] Wohlers | History of additive manufacturing[END_REF]. The principle of this process is similar to the FDM, while the feed material is in form of powder rather than filament and the applied material is usually metal. The powder feeding nozzle is coaxial with the laser beam, which melts the powder when it is propulsed and fed to the melt pool. This process is mainly used to repair complex and expensive parts instead of replacing them. It is suitable in the aeronautic field.

• Selective Laser Sintering (SLS) SLS was developed at the University of Texas at Austin at the mid-1980s [START_REF]Selective laser sintering[END_REF]. This is a powder-bed based process. It uses a laser as power source to bind the powder material, creating a solid structure. The used powder can be single-component, with only melting of the outer surface of particles. Two-component powders are mostly used, in which the powder with low melting point is melted and bind another non-melted powder to form solid parts. A wide range of materials can be used in SLS, including polymers and metals such as steel, titanium, alloy mixtures and composites. Ceramics are also used in this application.

• Selective Laser Melting (SLM) or Laser Beam melting (LBM) SLM started in 1995 at the Fraunhofer Institute ILT in Aachen [START_REF]Selective laser melting[END_REF]. This process is similar to SLS, but the powder is fully melted by an energy concentrated laser in SLM rather than sintered as in SLS. With SLM, homogeneous and fully dense parts can be obtained. Good mechanical properties are achieved after solidification. Another similar process is Electron Beam Melting (EBM), where an electron beam is used instead of a laser beam in order to melt the material.

The availability of different AM processes offers users the flexibility to choose suitable process according to the chosen material, construction type, desired precision and productivity. In the following, we will focus on SLM and present this process with more details.

Selective Laser Melting 1.2.1 Process principle

The principle of SLM is illustrated in Figure 1.2 with (a) an inside view of a real machine and (b) a schematic drawing of the process. The process is usually under the protection of a shielding gas, such as argon, in order to avoid oxidation. The emitted laser beam is focalized by lenses. A scanning mirror is used to deflect the laser beam toward the powder bed. By changing the deflection angle, the movement of laser spot in the XY plane can be controlled and high scanning speed can be achieved. Principal steps in SLM are described in Figure 1 • Initial step: Put the substrate at the bottom of build chamber, prepare process parameters and gas environment; (Note: see Figure 1.3 for steps 1 to 3)

• Step 1: The recoater arm (Figure 1.2) shoves the powder bed and deposits a layer of powder with controlled thickness ∆h p with respect to the top surface of substrate or previous consolidated layer;

• Step 2: The laser beam scans the powder according to a trajectory predefined by a 3D CAD model. The powder is melted and solidified to form a condensed layer with layer thickness H app , corresponding to a slice of the part (newly added red domain);

SELECTIVE LASER MELTING

• Step 3: The powder dispenser piston moves upward a distance ∆h up and the build piston moves downward a distance ∆h down . The powder of next layer is prepared to be deposited by the recoater arm;

• Repeat steps 1 to 3.

Here ∆h down is usually equal to the average height H app of constructed tracks in a layer when the process attains a stable regime. However, this height is smaller than the powder layer thickness ∆h p due to the shrinkage from powder to condensed parts. On the other hand, ∆h up is larger than ∆h down as the feeding powder should firstly fulfill the hollow caused by shrinkage. In addition, the powder may be compacted by the recoater arm when it is deposited on the previous layer in order to control the porosity. Additional attentions should be paid to ∆h down and ∆h up for the first several layers in order to obtain well attached tracks and better stabilize the track height. Two powder dispenser pistons may be used together in order to increase the deposition efficiency. In this case, the powder dispenser pistons move up alternatively in two successive layers and the recoater arm just needs to carry out a one-way movement rather than a round trip in each layer.

The procedure of layer construction is repeated to construct a 3D part layer by layer. The new molten and consolidated layer should be bound to the previous one to avoid mechanical defects. Besides the advantage of geometrical flexibility like any other additive manufacturing technologies, high relative density up to 99.9% can be obtained with SLM [START_REF] Yap | Review of selective laser melting: Materials and applications[END_REF], which means better mechanical properties. However, energy is so focused that process becomes more difficult to control. High thermal gradients usually exist and result in defects such as cracks. Thus, in order to obtain parts with high mechanical performance, the following process parameters should be carefully optimized:

• laser power, P L , nominal spot radius, r L and energy distribution;

• scanning speed, v L ;

• hatch spacing, ∆y (distance between two successive adjacent tracks within a given layer);

• powder layer thickness, ∆h p ;

• preheating temperature, T pre ;

• scanning and superposition strategies for 2D and 3D construction.

Laser, materials and powder

• Laser

The type of laser used in SLM is considered as an important factor affecting the properties of final parts. The choice of laser type is based on several considerations, including the power needed to melt the material, the spot size for the printing precision, the laser wavelength for the material absorption or the operating mode for a better heating. Two types of laser beams are mostly used in SLM, the Yb:YAG and CO 2 laser. The main difference between the Yb:YAG and CO 2 laser is the wavelength, to which attention CHAPTER 1. INTRODUCTION should be paid in application. The small wavelength of Yb:YAG laser (1.064 µm) makes it more suitable for metals than CO 2 laser (10.6 µm) as metals are highly reflective and CO 2 laser can not easily penetrate metals [START_REF] Yap | Review of selective laser melting: Materials and applications[END_REF]. On the other hand, Yb:YAG has difficulty to heat organic materials due to their weak absorption. For example, the light can simply pass through glass without heating it. However, the Yb:YAG lasers have higher beam quality and smaller spot size. Moreover, this type of laser is easier to manipulate and control. Another important point is the operating mode of laser, continuous or pulsed. Chou et al. [START_REF] Chou | Additive Manufacturing of Al-12Si Alloy Via Pulsed Selective Laser Melting[END_REF] show that SLM with pulsed laser applied to Al-12Si alloy allows a better control of heat input. Finer microstructure, higher density and hardness are obtained. In order to increase the production rate, some machines provide multiple lasers like the one with 4 lasers shown in Figure 1.2(a).

• Materials A wide range of materials is applicable to SLM. Studies focus principally on three types of metals: iron, titanium and nickel based alloys [START_REF] Yap | Review of selective laser melting: Materials and applications[END_REF] as shown in Table 1. Other aluminium, copper, magnesium, tungsten Table 1.1: Materials investigated in SLM, according to [START_REF] Yap | Review of selective laser melting: Materials and applications[END_REF] However, in spite of the attractive construction flexibility offered by SLM, few studies are conducted with ceramics due to the difficulties encountered in process, such as cracking and weak absorption of laser energy. The application of SLM to Al 2 O 3 -ZrO 2 ceramics is investigated in Fraunhofer Institute ILT by Wilkes et al. [START_REF] Wilkes | Rapid Manufacturing of Ceramic Parts by Selective Laser Melting[END_REF] and Hagedorn et al. [START_REF] Hagedorn | Net shaped high performance oxide ceramic parts by selective laser melting[END_REF]. By using a preheating laser, high strength part of almost 100% density without cracking was obtained. Nevertheless, poor surface quality was reported, caused by the large melt pool size due to the high preheating temperature.

• Powder

In SLM process, the initial state of the material is in the form of powder particles. The particle geometry results from powder processes. One of them is gas atomization, leading to a spherical geometry. The mechanical properties of final parts depend a lot on the properties of powder, including particle size distribution, morphology, porosity (or relative density), chemical composition, flowability and thermal properties [START_REF] Slotwinski | Characterization of Metal Powders Used for Additive Manufacturing[END_REF]. The porosity of powder bed is usually between 45% and 60%. It depends on the particle size One important issue in SLM is the powder recyclability in the point of economical and sustainable view. In order to ensure repeatable productions and desired mechanical properties, the characterization and chemical composition of recycled powder should not be largely affected during SLM. Strondl et al. [START_REF] Strondl | Characterization and Control of Powder Properties for Additive Manufacturing[END_REF] find that the new and recycled Ti-6Al-4V powder (Figure 1.4) in SLM has limited differences in processing and mechanical performances. The largest influence on the mechanical properties may be related to the modification of oxygen concentration during SLM. This study proves the feasibility of powder recycling in SLM.

Applications

Currently, SLM process is still more expensive than conventional processes such as casting and machining, due to the high price of powder and low production rate. Consequently, the application of SLM process is limited to parts with high value or requiring complex geometries.

The following are some typical applications of SLM.

• Aeronautic parts requiring high mechanical performance Aeronautic parts are usually of high value and with complex geometries, making it a suitable domain for the application of SLM. Jet engines have been realized by Monash University as shown in Figure 1.5(a). SLM is particularly attractive for parts working at high temperature like turbine blades and turbocharger rotors. These types of parts require good and reliable performance in creep, corrosion and thermal shock resistance. Nickel based superalloys like Inconel 718 are very suitable for this application. However, the excellent mechanical properties means also manufacturing difficulties. Machining of these parts often suffer from high wear rate of cutting tools, leading to high production price. The application of SLM to this type of superalloys can benefit the aeronautic industry as it provides more design freedom and keeps good mechanical properties. In addition, less or even no waste of expensive superalloys can be achieved. (c) 316L lattice structure for mechanical tests [START_REF] Alsalla | Fracture toughness and tensile strength of 316L stainless steel cellular lattice structures manufactured using the selective laser melting technique[END_REF]; (d) 316L resistojet with channels [START_REF] Romei | Manufacturing of a high-temperature resistojet heat exchanger by selective laser melting[END_REF].

• Medical orthopedic and dentistry

The application to medical orthopedic (Figure 1.5(b)) and dentistry is specially of interest due to the complexity and uniqueness of implants or dents. An interesting application of SLM is to make implants with controlled porosity. Li et al. [START_REF] Li | 316L Stainless Steel with Gradient Porosity Fabricated by Selective Laser Melting[END_REF] use SLM to fabricate pore gradient 316L stainless steel parts. They show the feasibility to produce implants with high porosity regions for the tissue growth while maintaining good mechanical properties regions with low porosity.

• Light weight structures Materials with high specific modulus (stiffness to weight ratio) and specific strength (strength to weight ratio) is always demanded in some applications like aerospace. This requirement drives the research on light weight structures. One possible way is to produce parts with controlled porosity while keeping equivalent stiffness and strength as bulk parts. The use of SLM process to produce lattice structures are investigated with different materials such as 316L (Figure 1.5(c) [START_REF] Alsalla | Fracture toughness and tensile strength of 316L stainless steel cellular lattice structures manufactured using the selective laser melting technique[END_REF]), Ti-6Al-4V [START_REF] Choy | Compressive properties of functionally graded lattice structures manufactured by selective laser melting[END_REF] and Al-Si10-Mg [START_REF] Maskery | A mechanical property evaluation of graded density Al-Si10-Mg lattice structures manufactured by selective laser melting[END_REF].

• Parts with cooling channels The design freedom with SLM process also benefits the fabrication of desired cooling channels, resulting in better control of cooling rate. This can be applied to injection moulding or aerospace like the resistojet heat exchanger shown in Figure 1.5(d).

PHYSICAL PHENOMENA IN SLM

Physical phenomena in SLM

Several physical phenomena are involved in SLM process, including heat and mass transfer, fluid dynamics in the melt pool, thermally induced stresses and microstructure evolution [START_REF] Markl | Multiscale Modeling of Powder Bed-Based Additive Manufacturing[END_REF].

The surface and mechanical properties of final parts are results of these phenomena. For instance, the heat transfer affects the temperature distribution, producing different stress distribution, which has direct influence on the formation of cracks and deformation. The track morphology and surface roughness are strongly related to the mass transfer and the fluid dynamics in the melt pool. The microstructure evolution is principally controlled by the heat transfer. The multiple physical phenomena in SLM result in difficulties of process understanding and control.

In addition, these phenomena occur at different length scales, varying from the periphery of the laser spot, to the melt pool, the solidified track and up to the whole part. This complicates the process and leads to the challenge of complete understanding of the effect of these phenomena. For example, the heat transfer around the laser spot is different to that in the melt pool and that in the whole part. On the other hand, the phenomena at small scale can have an influence up to the part scale. As an example, the mechanical properties of final part depend on the porosity and microstructure formed at the scale of the melt pool.

Consequently, a review on different phenomena in SLM is necessary both for process control and numerical modeling. These will be discussed in details in the following sections.

Heat and mass transfer

The driving force of the process is the input laser energy by irradiation. The heating and melting of powder is by the absorption of laser energy. The laser-powder interaction is very complicated due to the presence of shielding gas, the complex powder packing and sometimes the additional absorbers (see Section 5.1.1) for ceramics with low absorption. When a laser flux irradiates a powder bed, photons can be partially absorbed, reflected or transmitted by powder particles. Powder particles are heated by the absorption of photons. The reflections can be multiple between several particles, leading to a heating region larger than the nominal size of laser spot and also deeper than for bulk materials [START_REF] Gusarov | Modelling of radiation transfer in metallic powders at laser treatment[END_REF]. On the other hand, the transmission can be also important for certain materials transparent to laser like Al 2 O 3 and ZrO 2 (to Yb:YAG laser). This results into a much deeper penetration of the laser into the powder bed and, eventually, the consolidated substrate.

For a particle exposed to laser beam, the principal energy received is in the form of absorption. The absorbed energy is then redistributed by several heat transfer phenomena. The first is by conduction, through material or shielding gas. The former depends on the thermal conductivity of material while the latter is negligible due to the weak thermal conductivity of shielding gas. The conduction in material can be distinguished into intra and inter particles before melting. The intra-particle conduction can be considered as thermal conduction in bulk material. The inter-particle conduction is largely influenced by the sintering grade. Before sintering, particles are in point contact while surface contact occurs when sinter necks are formed, leading to a more efficient conduction [START_REF] Gusarov | Model of Radiation and Heat Transfer in Laser-Powder Interaction Zone at Selective Laser Melting[END_REF]. However, considering that the melt-CHAPTER 1. INTRODUCTION ing of powder is almost instantaneous, this transition is so rapid that such phenomena have probably no influence. The second phenomenon is by convection. In the powder bed, the convection firstly takes place in interstices fulfilled with shielding gas. For the same reason as previous, this can reasonably be neglected. Once the powder is melted, convection occurs principally in the melt pool. It is very important under a significant driving force like surface tension (including Marangoni effect), thus leading to the heat redistribution. The third phenomenon is by radiation, which has a biquadratic dependence on temperature. The last heat transfer phenomenon is the evaporation when the temperature reaches the material boiling point. This may occur when the laser power is too high or too concentrated, at low scanning speed or for a combination of these parameters. The evaporation can take away a lot of energy and cool the melt pool. The mass transfer is by the ejection of solid or liquid material in the form of powder denudation or liquid spattering, respectively. It can affect the track morphology and result into void structure. The denudation takes place around a track by clearing of powder particles, leading to a void adjacent area, as shown in Figure 1.6. By a detailed study with titanium and steel alloy powder, Matthews et al. [START_REF] Matthews | Denudation of metal powder layers in laser powder bed fusion processes[END_REF] find that the dominant driving force for denudation is the entrainment of particles by surrounding gas flow, due to the Bernoulli effect induced by the evaporation at the melt track center. An increase of denuded zone width with decreasing ambient pressure was observed. Both the powder denudation and spatter behavior in SLM are investigated by Liu et al. [START_REF] Liu | Investigation into spatter behavior during selective laser melting of AISI 316L stainless steel powder[END_REF] with 316L stainless steel. Their results show that the energy input affects the size, scattering and jetting height of spatter, usually in spherical 1.3. PHYSICAL PHENOMENA IN SLM shape. The spatters are furtherly classified into three types by Wang et al. [START_REF] Wang | Mechanisms and characteristics of spatter generation in SLM processing and its effect on the properties[END_REF] according to their morphologies. They concluded that these three types are mainly caused by recoil pressure due to the evaporation, Marangoni effect and heat effect in melt pool.

Melt pool dynamics

Melt pool dynamics is mainly driven by surface tension, Marangoni effect and recoil pressure. The influence of buoyancy force can be considered as minor due to the small melt pool size in SLM, even smaller than in welding [START_REF] Kruth | Consolidation phenomena in laser and powder-bed based layered manufacturing[END_REF]. It is crucial to understand the effects of these forces on the melt pool dynamics in order to obtain tracks with desired morphology.

The surface tension has tendency to smooth the surface of melt pool and has different effects on the track shape depending on the scanning speed. At low scanning speed, the melt pool is usually continuous. However, when the velocity is increased, the ratio of length and circumference of melt pool attains the limit of Plateau-Rayleigh instability [START_REF] Gusarov | Heat transfer modelling and stability analysis of selective laser melting[END_REF]. In order to minimize surface energy, the long melt pool is interrupted into fragments, which get spheroidized due to surface tension, as shown by Li et al. [START_REF] Li | Balling behavior of stainless steel and nickel powder during selective laser melting process[END_REF] in Figure 1.7. This balling effect is observed and analyzed by Gu and Shen [START_REF] Gu | Balling phenomena during direct laser sintering of multi-component Cu-based metal powder[END_REF][START_REF] Gu | Balling phenomena in direct laser sintering of stainless steel powder: Metallurgical mechanisms and control methods[END_REF]. Control methods are proposed for different kinds of balling effect [START_REF] Gu | Balling phenomena in direct laser sintering of stainless steel powder: Metallurgical mechanisms and control methods[END_REF].

The Marangoni effect is considered to have important consequences both on the melt pool geometry and the temperature inside the melt pool [START_REF] Chan | A two-dimensional transient model for convection in laser melted pool[END_REF][START_REF] Yuan | Molten pool behaviour and its physical mechanism during selective laser melting of TiC/AlSi10Mg nanocomposites: Simulation and experiments[END_REF]. It is caused by the surface tension gradient at the gas/liquid interface. The resulting convection flow is driven from low toward high surface tension zones. The surface tension depends on the temperature or the solutal concentration, thus a temperature or solutal gradient at the gas/liquid interface leads to the Marangoni effect and induces fluid flow.

As it has been mentioned in the previous section, the recoil pressure induced by evaporation is responsible of the powder denudation and liquid spattering. Moreover, due to this recoil pressure coupled with the Marangoni effect, additional liquid motion results into a deeper penetration of melt pool, leading to the keyhole formation [START_REF] Markl | Multiscale Modeling of Powder Bed-Based Additive Manufacturing[END_REF]. Gas can be trapped into the melt pool, which is solidified rapidly before gas escapes, thus leading to porosity in the final part. This phenomenon is detrimental as mechanical performance is decreased with the presence of porosity. Aboulkhair et al. [START_REF] Aboulkhair | Reducing porosity in AlSi10Mg parts processed by selective laser melting[END_REF] investigate different scanning strategies. They found that the keyhole induced porosity increases with the augmentation of scanning speed. The decrease of porosity is observed when they scanned the same layer twice.

Thermally induced stress

In SLM, laser heating is restricted to a small region around the laser spot, leading to a large thermal gradient and cooling rate. This induces considerable residual stresses and deformations in the final parts, leading to defects such as cracks and distortions. Elimination of cracks is a difficulty in SLM, especially for ceramics. Figure 1.8 shows the cracks observed with Al 2 O 3 -ZrO 2 system. It is formed (a) perpendicularly to the scanning direction and also (b) from the bottom of the melt pool as observed in the transversal cut view.

Two mechanisms of the residual stresses are analyzed by Mercelis and Kruth [START_REF] Mercelis | Residual stresses in selective laser sintering and selective laser melting[END_REF]. The first one is called the Temperature Gradient Mechanism (TGM), as illustrated in Figure 1.9.

It is caused by the large thermal gradient around the laser spot. Considering a thin sheet heated by a laser beam at the top, the neighbouring cold region will restrict the expansion of the heated region, thus inducing elastic compressive strains. Plastic strains are accumulated when it attains the yield strength. Then the compressed region starts shrinking when it is cooled and a bending angle towards the laser beam develops. The second one is related to shrinkage during the cooling of the melted region. The induced contraction is inhibited by the cold region and thus results into tensile stress in the heated region and compressive stress in the cold region.

Figure 1.9: TGM induced residual stresses [START_REF] Mercelis | Residual stresses in selective laser sintering and selective laser melting[END_REF] These two mechanisms indicate the approach to avoid defects by reducing the thermal gradient and the cooling rate. One solution is to use a second assisting laser for preheating. A CO 2 laser is employed by Hagedorn et al. [START_REF] Hagedorn | Net shaped high performance oxide ceramic parts by selective laser melting[END_REF], which covers the whole part with a preheating temperature of 1600 • C , for Al 2 O 3 -ZrO 2 eutectic ceramic. Cracks are prevented and high flexural strength is obtained.

Microstructure and texture

In SLM process, due to local melting and high scanning speed, rapid solidification can be achieved. This has a very important effect regarding the microstructural evolution, which induces large differences compared to conventional processes. It is well known that the SLM process can produce very fine and non-equilibrium structures [START_REF] Thijs | Fine-structured aluminium products with controllable texture by selective laser melting of pre-alloyed AlSi10Mg powder[END_REF]. Figure 1.10 shows microstructures of the Al 2 O 3 -ZrO 2 eutectic system, for which fine eutectic structure under micrometer is obtained. The different lamellar spacings at the center and at the border of the track are due to the evolution of the local solidification conditions. Similar to welding, solidification in SLM is directional, leading to morphological and crystallographic texture [START_REF] Thijs | Fine-structured aluminium products with controllable texture by selective laser melting of pre-alloyed AlSi10Mg powder[END_REF]. It is observed that the grain growth direction is perpendicular to the melt pool boundary (isotherm) and directed towards the center of the melt pool as shown in Figure 1.11. Thus grains are elongated in the direction of the thermal gradient, which is oriented to the building direction and largely influenced by the scanning speed. Consequently, the scanning speed, scanning and superposition strategies should be used in order to control the orientation and texture of grains.

Modeling of SLM process

The multiple physical phenomena in SLM are not sufficiently understood. This results in poor control of the process and increased efforts are needed to produce parts with desired geometrical and mechanical properties [START_REF] Markl | Multiscale Modeling of Powder Bed-Based Additive Manufacturing[END_REF]. The lack of process understanding is related to the challenges in the observation and measurement as several phenomena occur in a short duration and localized in a small region, such as the heat and mass transfer during the interaction between laser and material. Therefore, numerical modeling becomes an important approach to get a detailed insight of this process. In the following, different modeling scales and corresponding methods used in the literature will be briefly presented.

Micro modeling at track scale

Many physical phenomena can be investigated at track scale, such as heat transfer, melt pool dynamics, formation of solidification structures, stress field and so on. These phenomena are usually coupled together, leading to complicated modeling. Two approaches can be furtherly chosen, depending on whether powder particles are modeled [START_REF] Markl | Multiscale Modeling of Powder Bed-Based Additive Manufacturing[END_REF]. Modeling at this scale is usually limited to a few elementary tracks, due to intensive computing requirement, especially when considering approaches including the representation of powder particles.

• Powder approach

Modeling at powder scale provides a detailed study of physical phenomena in AM processes based on powder bed. The effect of laser-powder interaction, the formation of residual porosity and the development of track can be investigated with this approach. Mesh approaches are used in all modeling methods [START_REF] Markl | Multiscale Modeling of Powder Bed-Based Additive Manufacturing[END_REF]. A fine mesh at the scale of micrometers [START_REF] Körner | Mesoscopic simulation of selective beam melting processes[END_REF][START_REF] Khairallah | Mesoscopic simulation model of selective laser melting of stainless steel powder[END_REF] is necessary in order to provide a good representation of particle geometry.

The first issue raised up in this approach is the generation of powder bed. Körner et al. [START_REF] Körner | Mesoscopic simulation of selective beam melting processes[END_REF] and Shi and Zhang [START_REF] Shi | Simulation of random packing of spherical particles with different size distributions[END_REF] use the rain drop model [START_REF] Meakin | Restructuring effects in the rain model for random deposition[END_REF] to generate a 2D and 3D powder bed, respectively. Note that this approach is also available in Cemef, as applied in the thesis of Zouaghi [START_REF] Zouaghi | HIP of Stainless Steel 316L Considered at the Mesoscopic Scale: Numerical Modelling and Experimental Characterization[END_REF] to generate representative elementary volumes in the perspective of simulating hot isostatic pressure compaction applied to powder. As schematically shown in Figure 1.12(a), each particle is deposited separately. In analogy to a rain drop, the first contact of a falling particle is localized. Then it is rotated downwards until steady state. However, the relative density of powder bed (Figure 1.12(b)) generated by this method is approximately 74% in 2D and 60% in 3D, which is higher than the measured value between 45% and 60% (Figure 1.12(d)). Thus, particles are removed until desired porosity, causing non-physical holes in the powder bed (Figure 1.12(c)). Beside this, the difficulty of parallelization makes this method very time consuming for 3D cases. Ammer et al. [START_REF] Ammer | Simulating fast electron beam melting with a parallel thermal free surface lattice Boltzmann method[END_REF] and Yan et al. [START_REF] Yan | Multiphysics modeling of single/multiple-track defect mechanisms in electron beam selective melting[END_REF] use the discrete element method to generate powder bed. Each particle is assumed to move in a bounded space, under gravity, normal and tangential forces acting on each other. This method focuses on the relative density and the process stops once it is attained. The advantage of this method The melting of powder and melt pool dynamics result in significant surface evolution, thus requiring efficient surface tracking methods. Consequently, volume of fluid and level set methods are often used [START_REF] Markl | Multiscale Modeling of Powder Bed-Based Additive Manufacturing[END_REF]. The volume of fluid approach uses a so called fraction function for each cell. When this function is between 0 and 1, the cell is considered to be at the interface. A refinement process is required in order to track precisely the interface. The level set method uses a distance function and the interface can be directly represented by the zero iso-distance. It will be detailed later in this document as we use it in our approach.

• Continuum approach

The representation of powder particles requires a lot of mesh elements, leading to prohibitive computing charge. This is why the continuum approach using Finite Element (FE) method is usually employed. It considers the powder bed as a continuum. Under this assumption, most of the studies neglect the condensation of powder [START_REF] Gusarov | Heat transfer modelling and stability analysis of selective laser melting[END_REF][START_REF] Hodge | Implementation of a Thermomechanical Model for the Simulation of Selective Laser Melting[END_REF] and the resulting melt pool dynamics [START_REF] Yuan | Molten pool behaviour and its physical mechanism during selective laser melting of TiC/AlSi10Mg nanocomposites: Simulation and experiments[END_REF][START_REF] Li | A three-dimensional numerical model for a convection-diffusion phase change process during laser melting of ceramic materials[END_REF]. Hence, track morphology is not accessible, neither the residual porosity.

CHAPTER 1. INTRODUCTION

Considering the assumption of a continuous "powder material", its properties should take into account the presence of porosity in order to well simulate heat transfer. For the density and volumetric enthalpy, they can be simply averaged by volume fraction. However, for the thermal conductivity, special models or experimental calibrations should be used. A detailed review of effective thermal conductivity models is presented by Antwerpen et al. [START_REF] Van Antwerpen | A review of correlations to model the packing structure and effective thermal conductivity in packed beds of mono-sized spherical particles[END_REF], for powder beds packed by particles with the same size. The powder bed can be considered as a heterogeneous medium where discrete solid particles are packed randomly in a continuous gas phase. The effective conductivity of powder bed depends not only on the porosity, but also on the shape of particles, their sizes and their packing. It also results from convection exchanges through the gas and from radiation between powder particles. The low effective thermal conductivity of powder should be emphasized compared with bulk material.

Macro modeling at part scale

Modeling at the scale of part focuses generally on heat transfer and thermal induced deformation, usually based on Finite Element Method. This approach is particularly attractive to industry to design heat extraction in the presence of additional supports in order to control/master part deformation. Commercial software products are provided by several companies, such as Virfac of Geonx, Netfabb of Autodesk, exaSIM of 3DSIM, Simfact of MSC Software and the additive manufacturing package of ESI group. In order to reduce the computation charge, elements can be activated progressively according to the fusion of powder, leading to a simulation with "additively" increased active elements. This treatment is used by Chiumenti et al [START_REF] Chiumenti | Numerical modelling and experimental validation in Selective Laser Melting[END_REF] and simulations with large number of layers can be achieved. Similar to the continuum approach at micro modeling, the powder bed is treated as a continuous medium and continuous heat source models are often used, in surface or volume type. For example, the heat input is considered to be homogeneously distributed in a volume as modeled by Chiumenti et al. [START_REF] Chiumenti | Numerical modelling and experimental validation in Selective Laser Melting[END_REF] and Zhang et al. [START_REF] Zhang | Macroscopic thermal finite element modeling of additive metal manufacturing by selective laser melting process[END_REF]. In addition, in order to furtherly accelerate the computation, the equivalent energy deposition in track line or layer can be used [START_REF] Zhang | Macroscopic thermal finite element modeling of additive metal manufacturing by selective laser melting process[END_REF]. The representation of melt pool surface is not taken into account, thus no information about the solidified track surface is accessible. Under these approximations, the time step and minimum element size can be increased, giving an access to simulations with larger domain and longer duration. However, the element size is restricted by the layer thickness. As the present work focuses on the modeling at track scale, we will not discuss the macro modeling in more details hereafter. Only some results obtained at Cemef will be presented in order to give a panorama of the developed models at Cemef and the possible coupling of these two models in the future.

Challenges

The modeling of SLM is a challenging issue due to the complexity and interplay of multiphysic and multiscale effects. A complete modeling taking into account all phenomena at the scale of powder particles with simulation extended to the whole part is almost impossible. Just 1.5. CONCLUSIONS consider the fabrication of a 1 cm 3 cube with scanning speed v L = 1 m • s -1 , layer thickness H app = 50 µm and a hatch distance ∆y = 50 µm, the total scanning distance is 400 m and the scanning time is 400 s. Markl et al. [START_REF] Markl | Multiscale Modeling of Powder Bed-Based Additive Manufacturing[END_REF] point out that the time step of most models at powder scale is limited to the scale of nanoseconds. This makes the simulation very time consuming and unachievable for most of scientific computing. Consequently, simulations are usually restricted either at the micro or the macro scale, but not both. In addition, some physical phenomena are neglected in each type of modeling. For instance, the part distortion is inaccessible in the micro scale while the melt pool dynamic is usually not considered in the macro scale.

Conclusions

The history and principle of Additive Manufacturing are briefly introduced, together with its advantages and challenges currently encountered in the production. Different AM processes are compared according to their different principles and applications. The SLM process is then detailed, regarding the workflows of process, the applied laser, the materials and the applications.

The multiple physical phenomena occurring in SLM are presented. The heat transfer phenomena through conduction, convection, radiation and evaporation are detailed. The mass transfer caused by material ejection and liquid spattering is considered to have an important effect on the track morphology. Concerning the melt pool dynamics, the surface tension, Marangoni effect and recoil pressure have significant influence on the melt pool shape and thus track morphology. The formation of cracks is related to the high thermal gradient in SLM process and possible elimination method by preheating is presented in the literature. The fine structure due to rapid solidification and the texture linked to the directional solidification are also specific characterizations of SLM. At the end, modeling methods of SLM process at different scales are compared and the challenges of numerical modeling are emphasized.

Résumé

Après une étude bibliographique des méthodes de modélisation à différentes échelles dédiées à la fabrication additive basée sur un lit de poudre, la modélisation du procédé SLM à l'échelle du cordon, basée sur les méthodes éléments finis et méthode level set est présentée. Dans ce modèle, la forme des particules de poudre n'est pas représentée, la poudre est considérée comme un milieu continu afin de réduire la charge de calcul. Les différents modèles de source thermique sont présentés parmi lesquels celui basé sur la loi de Beer-Lambert est retenu. Ce modèle considère l'absorption locale du matériau. Il est adapté à la faible absorption des céramiques. L'implémentation de ce modèle est validée par un cas test avec absorption variable dans l'espace.

L'influence des propriétés du matériau et des paramètres du procédé sur la forme du bain liquide est étudiée à travers plusieurs cas tests. Le résultat justifie l'importance de l'absorption. De plus, la distribution de l'absorption dans différents milieux peut éventuellement modifier la distribution de la source thermique et donc la forme du bain liquide. Les paramètres du procédé ont également une grande influence sur le bain liquide. En conséquence, dans le procédé réel, on peut adapter les paramètres du procédé pour obtenir un bain liquide adapté avec les dimensions souhaitées, permettant d'optimiser les propriétés finales.

STATE OF THE ART

In this chapter, after a general introduction of modeling scales and methods, the literature study focuses on the modeling of heat transfer, including heat source models and results obtained in literature. Then the modeling method used in this work is presented, particularly the level set method and multiphase homogenization. The energy conservation and heat source model based on Beer-Lambert law are detailed, together with their numerical implementation. At the end, several cases are conducted to validate the heat source model and investigate the influence of material properties and process parameters on heat transfer.

State of the art

Due to the similarity between powder bed based AM processes, the literature study is not limited to modeling of SLM process, but also includes SLS and EBM. As we have mentioned in the first chapter, modeling can be divided into two scales: the track scale and the part scale. In this section, models for heat transfer corresponding to different scales are firstly presented based on a literature survey.

Heat source model

Given our context of numerical modeling applied to a continuous "powder material", the modeling of heat transfer requires an adequate heat source model. One simple heat source is the 2D Gaussian distribution (in the transverse section of the laser beam) multiplied by the absorptivity, which is used by Frewin and Stott [START_REF] Frewin | Finite Element Model of Pulsed Laser Welding[END_REF] in welding. The same surface heat source is also used in SLM [START_REF] Labudovic | A three dimensional model for direct laser metal powder deposition and rapid prototyping[END_REF][START_REF] Yin | Simulation of temperature distribution in single metallic powder layer for laser micro-sintering[END_REF][START_REF] Li | Thermal behavior during selective laser melting of commercially pure titanium powder: Numerical simulation and experimental study[END_REF]. However, considering the interaction between laser and porous powder, the heat input is more complicated than in welding processes. A volumetric heat source is more suitable due to the multiple reflection, especially for ceramics [START_REF] Li | Comparison of volumetric and surface heating sources in the modeling of laser melting of ceramic materials[END_REF]. The Goldak heat source model [START_REF] Goldak | A new finite element model for welding heat sources[END_REF] could also be used, which was originally developed for welding [START_REF] Desmaison | Modélisation Numérique d'un Procédé de Soudage Hybride Arc / Laser En Approche Level Set : Application Au Soudage Multi-Passes de Tôles d'acier de Forte Épaisseur[END_REF][START_REF] Chen | Three Dimensional Cellular Automaton-Finite Element (CAFE) Modeling for the Grain Structures Development in Gas Tungsten / Metal Arc Welding Processes[END_REF] in order account for fluid flow in the melt pool. However, this model relies on parameters deduced from a lot of measurements of melt pool dimensions and it is not explicitly based on the material properties. Zäh and Lutzmann [START_REF] Zäh | Modelling and simulation of electron beam melting[END_REF] use a volumetric heat source by adding an attenuation function in the propagation direction, in addition to the surface Gaussian distribution. However, the attenuation function is obtained by approximation of experimental results with a polynomial function, rather than based on physics. Consequently, this model may not be suitable for all materials and processes.

At the scale of powder particles, the laser/powder interaction can also be modeled in order to introduce a heat source. The ray tracing method has been proposed to track the trajectory of each light ray, then deduce the energy absorbed and scattered by powder particles. Zhou et al. [START_REF] Zhou | Numerical simulation of laser irradiation to a randomly packed bimodal powder bed[END_REF] use a two flux model (Figure 2.1) based on the balance of forward (red flux) and backward flux (blue flux) in an infinitesimal thickness dz:

   dφ + dz = -(α + β)φ + + βφ - -dφ - dz = -(α + β)φ -+ βφ + (2.1)
where α and β are the absorption and scattering coefficients, respectively. This model is CHAPTER 2. HEAT TRANSFER applied with a Monte Carlo method to examine the influence of particle surface emissivity and population ratio of larger and smaller particles in a bimodal random packing structure.

dz z φ + βφ + αφ + φ - βφ - αφ - (a
They find that such a bimodal packing structure with less large particles than small ones leads to higher radiative heat flux level. This model is also used by Khairallah et al. [START_REF] Khairallah | Laser powder-bed fusion additive manufacturing: Physics of complex melt flow and formation mechanisms of pores, spatter, and denudation zones[END_REF]. They find that heating begins at the particle surface and then diffuses inward with ray tracing, while a continuous volumetric heat source heats the whole particle instantly. In addition, the melting is not uniform with ray tracing as it occurs firstly on the surface and partial melting can be observed. They conclude that a heat source based on ray tracing is more realistic than a continuous volumetric source model.

Despite the precision of ray tracing method, it requires the representation of powder particles, thus only applicable to the modeling at particle scale. Gusarov and Kruth [START_REF] Gusarov | Modelling of radiation transfer in metallic powders at laser treatment[END_REF] developed a model based on the radiation transfer equation as follows:

Ω • ∇φ = -(α + β)φ + β 4π 4π φ Ω P Ω → Ω dΩ (2.

Heat transfer in AM processes

The modeling of heat transfer focuses on the prediction of the temperature field. A precise modeling should be coupled with phase change process. The temperature prediction during the laser melting of material can be achieved by the resolution of the non-steady convectiondiffusion heat transfer equation. At the micro scale, melt pool shape can be simulated and the influence of hydrodynamics on heat redistribution can also be investigated. On the other hand, macro modeling usually neglects the melt pool shape and takes the temperature distribution in the whole part as priority. The effect of additional support structures on heat extraction is also an important issue covered at this scale. Khairallah and Anderson [START_REF] Khairallah | Mesoscopic simulation model of selective laser melting of stainless steel powder[END_REF] present a complete 3D model with a random packing bed of 316L stainless steel , using a hybrid finite element and finite volume formulation on an unstructured grid. Note that the dimension of simulated system is limited at 1000 × 300 × 135 µm 3 . Their simulation couples the thermal diffusion and hydrodynamics, together with temperature dependent properties and surface tension. The results show the interesting islands of liquid regions as shown in Figure 2.2. They relate this phenomenon to the Plateau-Rayleigh instability, which creates peaks and troughs.

Hodge et al. [START_REF] Hodge | Implementation of a Thermomechanical Model for the Simulation of Selective Laser Melting[END_REF] present a simple micro model using the energy conservation equation taking account of the phase change by Stefan-Neumann equation. Powder particles are not represented and there is no coupling with hydrodynamics. The temperature field is predicted, as well as the shape of melt pool. As illustrated in Figure 2.3(a)-(d), when the laser passes a round trip on the powder non supported by the substrate (bulk region), two separate melt pools (in the overhang and bulk region) are obtained. They conclude that this phenomenon is consistent with experimental results. It is caused by the insulating behavior of the unconsolidated powder under the overhang region. The temperature histories of bulk and overhang node are shown in Figure 2.4. The first peak of each curve corresponds to the arrival of laser Li et al. [START_REF] Li | A three-dimensional numerical model for a convection-diffusion phase change process during laser melting of ceramic materials[END_REF] compare the influence of different heat transfer models on the melt pool shape in the simulation of laser melting of refractory ceramics. Three models are compared, pure heat transfer, with latent heat and with latent heat and fluid flow (driven by Marangoni effect). Their results in Figure 2.5 show that in this context of numerical modeling at the scale of melt pool, the complete model predicts better the melt pool shape and the pure heat transfer gives a high error. However, this model does not take into account the condensation of powder and no track height is predicted. For the modeling at part scale, we present briefly the model developed by Zhang et al. [START_REF] Zhang | Macroscopic thermal finite element modeling of additive metal manufacturing by selective laser melting process[END_REF] at Cemef. Complex part geometry is made possible in this model by the use of the level CHAPTER 2. HEAT TRANSFER set method. The construction path is totally process-oriented as the laser trajectory is taken from G-code files, which is written in numerical control programming language and mainly used in computer-aided manufacturing. One highlight of this work is the computing acceleration by the introduction of layer fractions as shown in Figure 2.6(a), which is equivalent to divide the scanning path into small successive subdomains. Each layer fraction is heated within a time equivalent to the exposure time of a material point to the laser spot. Then the fraction is cooled to complete the total scanning time of this fraction. By using a conform mesh adaptation strategy, simulations of the temperature evolution with a good representation of part geometry can be achieved (e.g. turbine in Figure 2.6(b)) under a reasonable computation time. Recently, it is coupled with thermomechanical model, providing access to the formulation of stresses and distortions during part construction. More details will be given in Chapter 4. Literature results show the importance of hydrodynamics on heat transfer. In addition, the insulation effect of powder is a specialty of SLM. In the following, the modeling of heat transfer with a particular heat source model for ceramics will be introduced. It is coupled with hydrodynamics but coupled simulations will be demonstrated in the next Chapter when the melt pool dynamics is modeled.

MODELING METHOD

Modeling method

As it is difficult to couple the macro and micro scale modeling in a single model, the modeling scale should be chosen considering the physical phenomena of interest. In the present studies, the temperature field, the melt pool dynamics, the track shape and the potential crack formation meet the interest of the partner of C ÉFAL É project. The modeling introduced hereafter is at the scale of track formation, without the representation of powder particles. The modeling considers a whole system consisting of substrate and possibly the previously deposited layers, powder bed and gas. The FE method will be used, together with the Level Set (LS) method to capture the track surface. Principal assumptions should be mentioned before:

• the powder bed (particles with gas) is assimilated to a continuous medium;

• no powder projection (i.e. no powder loss) or evaporation will be taken into account;

• no residual porosity once the powder is fully melted (i.e, no release of dissolved gas to form pores; no gas entrapment). The present approach is based on preliminary activity on numerical simulation of welding hybrid arc-laser welding developed at Cemef by Hamide [2] and Desmaison [START_REF] Desmaison | Modélisation Numérique d'un Procédé de Soudage Hybride Arc / Laser En Approche Level Set : Application Au Soudage Multi-Passes de Tôles d'acier de Forte Épaisseur[END_REF][START_REF] Desmaison | A level set approach for the simulation of the multipass hybrid laser/GMA welding process[END_REF]. The system consists of two domains, named material (D 1 ) and gas (D 2 ), as shown in Figure 2.7.

The porous powder layer is deposited on the fully dense substrate. The densification by the melting of powder is modeled through the variation of apparent density (mass occupied by an unit volume including dense material and porosity). The shape of the melt pool and the final track is affected by the forces acting on the liquid, principally the surface tension, the Marangoni effect. The evaporation pressure is not considered in the present modeling. Note that this phenomenon may be not as important as processing metals due to the high boiling temperature (e.g. 3240 K for Al 2 O 3 [START_REF] Langstaff | Aerodynamic levitator furnace for measuring thermophysical properties of refractory liquids[END_REF]) of ceramics. The evolution of gas/material interface is tracked by a level set and mesh adaptation. All conservation equations are based on the CHAPTER 2. HEAT TRANSFER whole system with different material properties for gas and material domains as detailed hereafter.

Level set approach

The level set method is firstly introduced by Osher and Sethian [START_REF] Osher | Front Propagating with Curvature Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations[END_REF]. This method provides the access to the tracking of evolving object surface without parameterization of the object.

In addition, this method can easily follow the topological change of objects, e.g. splitting of a bubble. This property makes it largely used in computational fluid dynamics [START_REF] Sethian | Level set methods for fluid interfaces[END_REF]. Considering the surface evolution of the melt pool in SLM, this method is suitable in such a continuous model.

Taking a system composed of two domains, said gas (D 1 ) and material (D 2 ) in the context of SLM. The idea of LS method is to defined a signed distance function ψ with respect to the interface Γ (ψ = 0) of D 1 and D 2 (Figure 2.7):

ψ(x) =          -d(x) if x ∈ D 1 0 if x ∈ Γ d(x) if x ∈ D 2 (2.5)
where d(x) is the geometrical distance of point x to the interface Γ. The LS function should respect the Eikonal equation ∇ψ = 1. Consequently, the geometrical distance d(x) between two iso-contours ψ 1 and ψ 2 is equal to |ψ 1 -ψ 2 |. This condition needs to be systematically valid even during the update of ψ due to interface evolution. The normal direction of isosurface, in particular the one corresponding to ψ = 0, is then defined as:

n = ∇ψ ∇ψ (2.6)
As mentioned above, the level set ψ = 0 follows the mobile gas/material interface driven by a velocity field u. The transport equation is thus solved to update the level set function:

dψ dt = ∂ψ ∂ t + u • ∇ψ = 0 (2.7)
This solution gives access to the new position of the gas/material interface (ψ = 0). However, the updated field ψ does not necessarily respect the Eikonal equation. Therefore, ψ is usually reinitiallized after or during the transport. Two principal methods are proposed in the literature. The first one is to solve the convected reinitialization equation [START_REF] Ville | Convected level set method for the numerical simulation of fluid buckling[END_REF]:

∂ψ ∂ t + (u + bU ) • ∇ψ = bs (2.8)
where b is a coefficient usually chosen equal to the ratio between the mesh size and the time 2.2. MODELING METHOD step. s is the sign function of ψ and U is computed as:

U = s ∇ψ ∇ψ (2.9)
Note that the transport and reinitiallization of ψ is done in a single resolution in Eq.2.8. However, the convection of ψ in the whole domain may lead to numerical instabilities. In fact, only the interface position ψ = 0 is of interest. Thus a filtered distance function can be used while keeping the exact position of ψ = 0. The resolution of the convected reinitiallization equation provides a good interface representation. Nevertheless, the use of filter function does not give a real distance value, leading to the loss of precision in the computation of the normal direction or curvature.

Another method is more direct based on the geometrical consideration [START_REF] Shakoor | An efficient and parallel level set reinitialization method -Application to micromechanics and microstructural evolutions[END_REF]. After solving Eq.2.7, the iso-surface ψ = 0 can be obtained. Then the distance of each mesh node to this iso-surface can be computed geometrically. This method, firstly experimented at Cemef by Desmaison in his thesis [START_REF] Desmaison | Modélisation Numérique d'un Procédé de Soudage Hybride Arc / Laser En Approche Level Set : Application Au Soudage Multi-Passes de Tôles d'acier de Forte Épaisseur[END_REF], was later implemented by Shakoor et al. [START_REF] Shakoor | An efficient and parallel level set reinitialization method -Application to micromechanics and microstructural evolutions[END_REF] in Cimlib.

Homogenization of properties

The computation of global property can be carried out at two levels, including the Heaviside average and the multiphase homogenization as detailed hereafter.

Property averaged by Heaviside function

Material properties, such as density and conductivity, can have a difference of several order of magnitude between the two domains D 1 and D 2 . A sharp change of material properties at the interface ψ = 0 may lead to numerical instabilities. A transition zone is thus introduced around the interface ψ = 0 with a half thickness , where a Heaviside function H evolves continuously from 0 to 1:

H(ψ) =            0 if ψ < -(material domain D 1 ) 1 2 1 + ψ + 1 π sin πψ if | ψ | 1 if ψ > (gas domain D 2 ) (2.10)
The derivation of Heaviside function leads to the Dirac function δ:

δ(ψ) =            0 if ψ < -(material domain D 1 ) 1 2 1 + cos πψ if | ψ | 0 if ψ > (gas domain D 2 ) (2.11)
The global properties {χ} are averaged between two domains by the Heaviside function [START_REF] Chen | Three Dimensional Cellular Automaton-Finite Element (CAFE) Modeling for the Grain Structures Development in Gas Tungsten / Metal Arc Welding Processes[END_REF][START_REF] Desmaison | A level set approach for the simulation of the multipass hybrid laser/GMA welding process[END_REF], in particular in the transition zone around the D 1 /D 2 interface:

{χ} = H χ D 2 + (1 -H) χ D 1
(2.12) ; Figure 2.8: Heaviside and Dirac function here χ D i (i = 1 for material and i = 2 for gas) denotes the average properties in each domain homogenized between phases by their volume fractions. More details will be given in the next section. This treatment is used for density, enthalpy, thermal conductivity and viscosity. It aims at defining a set of equations averaged over the whole system.

Multiphase homogenization

The material domain (D 1 ) is a multiphase domain. A general description is introduced by dividing this domain into zone, structure and phase, as shown in Table 2.1, for the case of pure Al 2 O 3 . This approach provides a multi-level description of the material domain. For example, the zone notation can separate the powder bed (Z 1 ) and dense matter (Z 2 , including the substrate and consolidated tracks). Inside the powder bed, we can find the solid powder particles (S 1 ) and the porosity (S 2 ). At a smaller scale, the solid structure consists of the phase α-Al 2 O 3 . The advantage of this description is that the volume fraction of each phase, structure and zone can be tracked and the material properties of this domain can be homogenized based on different levels. Note that in this case, there is only one phase in each structure and the distinction between structure and phase may be not necessary. However, this description can be easily extended to other materials like the eutectic of Al 2 O 3 -ZrO 2 , where the phases α-Al 2 O 3 , monoclinic ZrO 2 and tetragonal ZrO 2 can exist in the solid structure (or eutectic structure) S 1 . The tabulation is simplified in the gas domain, as there is only one zone, structure and phase inside. The material properties in the domain can be averaged in the level of phase, structure or zone. For example, density and volumetric enthalpy are averaged between phases by their volume fractions, according to the concept of Representative Element Volume (REV) [START_REF] Rappaz | Numerical Modeling in Materials Science and Engineering[END_REF]: where χ β i is the intrinsic property of phase β i , g β i S j the volume fraction of phase β i in the structure S j , g S j Z k the volume fraction of structure S j in the zone Z k and g Z k D 1 the volume fraction of zone Z k in the domain D 1 . Note that for zone, structure and phase, the sum of their fractions is equal to 1:

χ D 1 = β i ∈D 1 g β i D 1 χ β i = Z k ∈D 1 S j ∈Z k β i ∈S j g Z k D 1 g S j Z k g β i S j χ β i (2.
Z k ∈D 1 g Z k D 1 = 1, S j ∈Z k g S j Z k = 1, ∀Z k , β i ∈S j g β i S j = 1, ∀S j (2.14)
Different average laws in different level can be applied to compute χ D 1 , based on physical considerations. For example, the thermal conductivity of material may be averaged between powder (Z 1 ) and dense matter (Z 2 ) in an arithmetic way, while special non-linear models can be used to compute the effective conductivity of powder.

The evolution of phases, structures and zones depends on the temperature. Their volume fractions can be tabulated with respect to the temperature, according to the melting and solidification path. This will be detailed in Section 2.4.1.

Heat source model

In the context of powder bed based AM process, as laser interacts with different media like powder and dense matter, a heat source model taking into account material properties is more suitable. Essentially, the local absorption coefficient has to be considered. In the following, we will try to model the heat source based on the Beer-Lambert law.

Material absorption

As presented by Hagedorn [START_REF] Hagedorn | Additive Manufacturing of High Performance Oxide Ceramics via Selective Laser Melting[END_REF], for an absorbing medium, its intrinsic optical index n can be expressed in complex form:

n = ϑ + iκ (2.15)
Here ϑ is the refractive index, which represents how much the light is refracted when entering into a material. This parameter is the ratio between the light speed in vacuum and in the material. κ is the extinction coefficient, which indicates the energy attenuation of an electromagnetic wave propagating in this medium, i. e. the energy absorbed by the material.

Considering an initial electromagnetic wave E 0 (z, 0), when it propagates (in z direction) in a CHAPTER 2. HEAT TRANSFER homogeneous bulk material with constant optical index, it can be written by: E(z, t) = E 0 (z, 0)e i(nωz/c-ωt) = E 0 (z, 0)e -κωz/c e i(ϑωz/c-ωt) (2.16) where ω is the angular frequency and c is the light speed. Note that the real part of E describes the decrease in propagation direction z, while the imaginary part corresponds to the sinusoidal oscillation. The optical intensity of a light wave φ is proportional to EE * (E * is the conjugate of E). Thus the Beer-Lambert law with constant α can be obtained as:

φ = C • E 2 0 e -2κωz/c = φ 0 e -αz (2.17)
where C is the proportion coefficient, φ 0 is the initial intensity and the absorption coefficient α, is defined as:

α = 2κω c = 4πκ λ L (2.18)
where λ L is the laser wavelength. Note that with a given laser of wavelength λ L , the absorption of material is directly related to its extinction coefficient κ. The latter depends both on the laser wavelength and the temperature [START_REF] Faure | Etude de l'interaction Rayonnement-Matière Dans Un Milieu Granulaire En Vue de l'application Au Procédé de Frittage Laser[END_REF].

Yb:YAG laser . This means that when a flux penetrates an absorption length of about 0.1 m (1/α), its intensity only decreases by 36.8%. This is a very low absorption compared with metals, which absorb laser radiation almost at the surface. However, the presence of impurities can largely modify the absorption coefficient. Lawrence [START_REF] Lawrence | An analysis of the beam interaction characteristics of selected lasers with an alpha-alumina bioceramic[END_REF] reports 1/α = 106 ± 7 µm for 99.4% pure alumina with 0.6% impurities of CaO, Fe 2 O 3 , K 2 O, MgO, SiO 2 and Ti 2 O. This corresponds to α = 8.85 ∼ 10.1 mm -1 , which 2.3. HEAT SOURCE MODEL is about 1000 times higher than the single crystal in [START_REF] Faure | Etude de l'interaction Rayonnement-Matière Dans Un Milieu Granulaire En Vue de l'application Au Procédé de Frittage Laser[END_REF]. In conclusion, it is very difficult to get a reliable value of the absorption coefficient as it depends a lot on impurity composition. In the process, absorbers [START_REF] Defillon | Modèle du bain de fusion de SLM sur céramique avec disparition de l'agent absorbant restreinte au bain[END_REF] like carbons are usually added in order to increase the heat efficiency and homogenize the absorption.

Volumetric heat source model

The absorption coefficient α can vary depending on the material state and the local phase fractions (liquid phase and/or different solid phases). The Beer-Lambert law in the form of partial equation applied to heterogeneous material (Eq.2.3) can be solved by combining the boundary condition φ = φ 0 at z = 0, leading to:

φ(r, z) = φ 0 (r, 0) exp - z 0 α(r, z)dl (2.19)
for a local position (r, z) considering a propagation of laser flux only in the z direction (Figure 2.7).

Assuming that the attenuation of laser energy in a thickness dz is totally absorbed by the material, we can then deduce a volumetric heat source model, which is equal to the flux variation in propagation direction z:

qL = - dφ dz = αφ 0 exp - z 0 αdl = αφ (2.20)
where the initial surface flux φ 0 is assumed to follow a radial Gaussian distribution with a standard deviation equal to half of the interaction radius r int :

φ 0 = φ(r, 0) = 2P L πr 2 int exp - 2r 2 r 2 int (2.21)
Here the interaction radius r int is used rather than the nominal laser radius r L due to the multiple reflection of laser flux in powder bed, leading to r int > r L [START_REF] Gusarov | Modelling of radiation transfer in metallic powders at laser treatment[END_REF]. This initial flux will be reflected at the material surface ψ = 0 with reflection coefficient R. Note that the surface reflection can also vary depending on the material state, as dense solid, liquid and porous powder have different reflection coefficients. However, an average effective coefficient will be considered in the following.

The volumetric heat source model can be finally expressed by:

qL (r, z) = (1 -R) • 2P L πr 2 int exp - 2r 2 r 2 int • α • exp - z 0 αdl (2.22)
At the Right Hand Side (RHS) of this equation, the four terms represent the reflection, the Gaussian distribution, the local distribution and the laser intensity attenuation during propagation, respectively. The local variation of absorption coefficient is taken into account in the integration along the z-direction in the last term. In the numerical implementation, the heat source deposition is truncated in the region r ∈ [0, 1.5r int ]. Consequently, the total power deposited in this region is 1 -e -9/2 = 98.9%P L (1 -R). Thus normalization is conducted to deliver the nominal net power P L (1 -R). The heat source is proportional to P L and R, while the influence of r int and α is more complicated. Hence iso-contours of qL under different r int and α are compared in Figure 2.10. The increase of r int always leads to shallower iso-contours, while it becomes narrower (e.g. iso-contour 0.4) near the center but wider (e.g. iso-contour 0.1) far away from the center. On the other hand, the decrease of α always results in narrower iso-contours, while its influence on the iso-contour depth is not monotonous. With smaller values for α coefficient, the region with high energy becomes shallower (e.g. iso-contour 0.4) and the one with low energy becomes deeper (e.g. iso-contours 0.1 and 0.2).

Energy conservation

The energy received by powder results in the increase of temperature and the phase transformation from solid to liquid. At the rear of melt pool, the heat extraction by conduction with surrounding solid cools and solidifies the melt pool, thus forming the track. The energy conservation solution should be coupled with the melting and solidification path by the evolution of thermodynamic variables related to each phase like fraction, density and specific enthalpy [START_REF] Saad | Temperature-based energy solver coupled with tabulated thermodynamic properties -Application to the prediction of macrosegregation in multicomponent alloys[END_REF]. Hence, after the resolution of energy conservation, we can obtain not only the temperature, but also these variables. This provides a detailed description of material state and other material properties like conductivity can be furtherly computed. In the following, we will firstly define melting and solidification path and then present the energy conservation equation coupled with it.

Phase transformation

When pure alumina begins to melt or solidify, the enthalpy increases or decreases respectively while keeping at a constant temperature -the melting point T m = 2054 • C. Phase transformation takes place at this temperature with a large enthalpy variation due to the release of latent heat. When proceeding to numerical modeling, this causes difficulties as one temperature can correspond to different enthalpy states. In order to avoid this problem, an artificial solid-liquid transformation interval is assumed between 2004 • C (∼solidus) and 2104 • C (∼liquidus). The volume fraction of solid phase in dense matter (g s Z 2 = g s S 1 g S 1 Z 2 ) evolves linearly from 1 to 0 with respect to temperature from 2004 • C to 2104 • C and inversely for the liquid phase (g l Z 2 = g l S 3 g S 3 Z 2 ). In the tabulation, g s S 1 and g l S 1 are always equal to 1, while g S 1 Z 2 (or g S 3 Z 2 ) changes from 1 to 0 (or 0 to 1). The melting of powder induces simultaneously the one-way condensation from powder to dense matter as interstitial porosity is fulfilled by the liquid. This leads to a significant shrinkage and thus the displacement of gas/material interface. Note that the heating rate in SLM is very high and this happens in a small time interval on the order of 100 µs [START_REF] Khairallah | Mesoscopic simulation model of selective laser melting of stainless steel powder[END_REF]. The resulted high convection velocity requires the decrease of time step according the law of Courant-Friedrichs-Lewy [START_REF] Markl | Multiscale Modeling of Powder Bed-Based Additive Manufacturing[END_REF]. Hence, we propose an enlargement of the condensation temperature range from 1804 • C to 2104 • C, in which the powder fraction (g Z 1 D 1 ) evolves CHAPTER 2. HEAT TRANSFER linearly from 1 to 0 and inversely for dense matter (g Z 2 D 1 ). This means a condensation step before melting. Considering the short duration of this step, the impact of this assumption on the prediction final bead shape should be quite limited as the bead shape essentially depends on fluid flow.

Figure 2.11 shows an example of zone evolution and a schema with corresponding phases, structures and zones in an Representative Elementary Volume (REV) [START_REF] Dantzig | Solidification[END_REF] in each state. In state 1 below the condensation temperature 1804 • C, there is only powder (Z 1 ). Heating until state 2 leads to the shrinkage and powder is partially condensed to dense matter (Z 2 ). When the temperature is higher than the solidus in state 3, liquid is formed. In this state, material can be directly cooled to state 6 (partial melting) or furtherly melted to attend state 4 (full melting) and then cooled to state 5. The dashed lines in (a) indicates the one-way transformation, meaning that the powder fraction can only be decreased while the fraction of dense matter can only be increased.

Conservation equation

The modeling of heat transfer in SLM consists of the resolution of energy conservation equation. This resolution provides the temperature distribution in space and its evolution in time. This step is crucial for the prediction of phase fractions, and furtherly the induced fluid dynamics in the melt pool and the stress distribution. The conservation equation can be established in the whole domain including gas and material:

∂{ρh} ∂ t + ∇ • {ρh}u -∇ • {λ}∇T = Q = qL -qr {ρh} = f (T ) (2.23)
where ρ is the density, h the specific enthalpy, λ the thermal conductivity. All the quantities are considered in their average formulation, noted {•}. Q is the net input power, including the heat source qL and the heat loss qr by radiation at the immersed gas/material interface Γ:

qr = δεσ r T 4 -T 4 ext (2.24)
where σ r = 5.67

× 10 -8 W • m -2 • K -4
is the Stefan-Boltzmann constant. ε is the emissivity and T ext the environment temperature. The surface loss is transformed into volumetric loss by the method of Continuous Surface Force [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF]. Boundary conditions can be of several types:

• Dirichlet condition (imposed temperature):

T = T imp on ∂Ω T (2.25) • Neumann condition (imposed flux): λ∇T • n Ω = q imp on ∂Ω q (2.26)
• Convection and radiation:

2.5. NUMERICAL RESOLUTION λ∇T • n Ω = h c + εσ r (T + T ext ) T 2 + T 2 ext (T -T ext ) on ∂Ω cr (2.27)
The second term in the Eq.2.23 can also be expanded to:

∇ • {ρh}u = u • ∇{ρh} + {ρh}∇ • u (2.28)
The first term at the RHS is related to the convection in Eulerian formulation. The second term vanishes if the material is incompressible (∇ • u = 0). However, in our modeling of SLM, the shrinkage from powder to dense matter, which is due to the increase of apparent density, leads to a non null value of ∇ • u. This will be discussed in the mass conservation in the next chapter.

An important issue is the relationship between T and {ρh}. Volume fractions and specific enthalpy of each phase can be obtained at given temperature T and computed by the temperature-enthalpy conversion where melting and solidification paths are predefined. Then with Eq.2.12 and 2.13, the global volumetric enthalpy {ρh} can be obtained. The non-linear effect of Eq.2.23 raises up due no the non-linear relationship between T and {ρh}, especially with the latent heat during solid-liquid transformation. The numerical resolution of this non-linear equation will be detailed in Section 2.5.2.

Numerical resolution

Heat source implementation

Parallelization strategy

As the volumetric heat source model in Eq.2.22 takes into account the variation of absorption coefficient α, the last attenuation term requires numerical integration. In order to integrate along one direction, a series of successive sampling points is required. However, this set of points is not provided by the triangle or tetrahedron mesh used in FE method. Thus, the integrand field (here absorption coefficient α) is interpolated to a virtual grid used only for the integration operation. After the integration on grid, the resulting integral field is interpolated back to the FE mesh. The schema in Figure 2.12(a) illustrates the relationship between the grid and the FE mesh for a 2D case. In order to interpolate from the FE mesh to the grid, for a given grid point x g i (i ∈ [START_REF] Hagedorn | Net shaped high performance oxide ceramic parts by selective laser melting[END_REF]NbGrid] where NbGrid is the number of grid points), we should firstly find the FE element E k (k ∈ [START_REF] Hagedorn | Net shaped high performance oxide ceramic parts by selective laser melting[END_REF]NbE] where NbE is the number of FE elements) containing this point Figure 2.12(a). This step is time consuming and adequate parallelization is required.

Cimlib is already parallelized with respect to the FE mesh as shown in Figure 2.12(a). If we keep this parallelization strategy, for each grid point, we should scan the elements in each partition to find the one containing it. This may be extremely slow if the distribution of grid points is not equivalent in each region occupied by FE mesh. For example, if the grid is totally covered by one partition, this partition will be highly charged while others have nothing to do. This parallelization strategy is consequently inefficient. A more efficient algorithm can be achieved just by inverting the parallelization strategy. Firstly we can assemble the FE mesh in each partition (red, blue and green FE elements in Figure 2.12(a)) to obtain a global mesh (solid black in Figure 2.12(b)) and broadcast it to each core, similarly for the integrand field α. Then the grid points can be partitioned equivalently in each core (red, blue and green grid points in Figure 2.12(b)). Thus, points in each partition can localize their containing elements independently, as every partition has the information of global mesh. This localization strategy will be detailed in the Section 2.5.1.3.

By contrast, the interpolation from the grid to the FE mesh is very quick. Due to the regularity of the grid, the localization of rectangle (2D) or cuboid (3D) containing the mesh node x j (Figure 2.12.(b)) can be simply carried out with the coordinates of x j , making it much faster than the search of finite element.

Interpolation and integration

The interpolation from FE mesh to grid can use the classical treatment in FE method [START_REF] Rappaz | Numerical Modeling in Materials Science and Engineering[END_REF]. The interpolation function N j x g i of each node j of element E can be evaluated at the grid point x g i as illustrated in Figure 2.13(a). The value of integrand field at this point is:

α x g i = D j=1 N j x g i α x j (2.29)
where D is the number of nodes of a FE element (3 for 2D and 4 for 3D). Several numerical methods exist for the 1D integration, such as trapezoidal rule and Simpson's rule. Although the Simpson's rule (∆ 5 x) is two orders more accurate than the S EF GH from grid to FE mesh. trapezoidal rule, the latter is used in the following treatment as it is more simple with sufficient precision. Once the integrand field α x g i associated to the grid is obtained, the numerical integration by trapezoidal rule is applied:

A x g i = x g i x g 1 α(x g ) dl = i-1 k=1 α x g k+1 + α x g k 2 ∆l (2.30)
where x g k and x g k+1 are successive grid points along the integration direction as illustrated in Figure 2.13(b). Finally, the integral field obtained on the grid is returned to the FE mesh by bilinear (2D) or trilinear (3D) interpolation:

A x j = 2 d i=1 N g i x j A x g i (2.31)
where d is the dimension of space and N g i x j the interpolation function (from grid to FE mesh) of grid point i evaluated at FE node x j . It can be computed as shown in Figure 2.13(c).

Resolution schema

The resolution schema of heat source distribution qL is shown in Figure 2.14. The localization (steps in cyan rectangle) of grid point (x g i ) in the containing element (E k ) is optimized in several levels. Once the mesh is updated or the region impacted by laser gets out of the grid, relocalization should be conducted. In this situation, for the first point (i = 1), the distance from this point to the FE element center is computed. Then the search of containing element E k begins with the nearest element E l . This process checks E l and its neighbors (not only the first patch) until find the containing element. For other grid points, the search is processed CHAPTER 2. HEAT TRANSFER successively. Each point x g i will check if the containing element E l of its neighboring grid

Start

Assemble mesh, α, create grid Loop over

x g i , i ∈ [1, NbGrid] Relocalization i = 1 Find E l near- est to x g i With E l , search E k , x g i ∈ E k Interpolation FE to grid Find E k in memory, x g i ∈ E k Find neigh- bor x g j of x g i x g j ∈ E l obtained i =
, i, j ∈ [1, NbGrid], k, l ∈ [1, NbE].
point x g j is already obtained precedently. If it is true, this element E l can be used as start value to search the element containing x g i . If not, the same process used for the first grid point will be conducted. Due to these treatments, the search efficiency can be highly increased.

On the other hand, if the relocalization is not required, the localization relationship is the same as the previous, which is saved in memory. Hence the containing element E k can be directly obtained. In this case, the computation time is largely decreased as the search of containing element is omitted.

FE discretization

The energy conservation equation in Eq.2.23 can be solved by Finite Element (FE) method. It is firstly transformed into a weak form by a variational formulation. This is carried out by multiplying Eq.2.23 by a test function w belonging to the Sobolev space H 1 (Ω) and integrate in the whole domain Ω. The second step is the discretization of Ω to Ω E by using finite element. In the following, simple linear element (P 1) will be used, which is a triangle in 2D and tetrahedra in 3D, with number of nodes D=3 and 4 in 2 and 3 dimensions, respectively. The problem can be finally expressed by:

Find T ∈ T = T : T ∈ H 1 Ω E , T = T imp on ∂Ω E T , so that for ∀w ∈ W = w : w ∈ H 1 Ω E , w = 0 at ∂Ω E T , satisfying : Ω E w ∂H ∂ t dV + Ω E w∇ • {ρhu} dV - Ω E w∇ • {λ∇T } dV - Ω E w Q dV = 0 (2.32)
where {ρh} is noted as H. Here the test function w of classical Galerkin type is used for the simplicity, although the implementation in the Cimlib is stabilized by SUPG method [START_REF] Hachem | Stabilized Finite Element Method for Heat Transfer and Turbulent Fows inside Industrial Furnaces[END_REF].

Eq.2.32 is the global system to be determined. In fact, it can be assembled by local matrix contributed by each element E at time t [START_REF] Saad | Numerical Modelling of Macrosegregation Formed during Solidification with Shrinkage Using a Level Set Approach[END_REF]:

R E i t = M E ij H t j -H t-∆t j + A E ij H t j + K1 E ij + K2 E ij T t j -F E i -Q E i = 0 (2.33)
where (i, j ∈ [1, D]) and the local matrices related to element E is (only convection boundary condition is considered):

M E ij = E 1 ∆t N i N j dV (2.34
)

A E ij = E N i u • ∇N j + N i N j ∇ • u dV (2.35
)

K1 E ij = E {λ}∇N i • ∇N j dV (2.36) K2 E ij = ∂E h c N i N j dV (2.37
)

F E i = ∂E h c T ext N i dS (2.38
)

Q E i = E QN i dS (2.39)
The local matrix is computed by looping all the volumetric and surface elements. It is assembled into a global matrix, leading to a final linear system to be solved. This system can be expressed as the same form as Eq.2.33 while it is applied to the whole domain rather than the element E and the global numbering will be used for node index i, j.

(R i ) t = M ij H t j -H t-∆t j + A ij H t j + K1 ij + K2 ij T t j -F i -Q i = 0 (2.40) CHAPTER 2.

HEAT TRANSFER

As we have mentioned before, this is a non-linear system due to the non-linear relationship between temperature and volumetric enthalpy. One classical resolution method of non-linear system is by Newton-Raphson iteration [START_REF] Rappaz | Numerical Modeling in Materials Science and Engineering[END_REF]. Considering the following system to be solved:

R(V ) = 0 (2.41)
Suppose that at iteration (ν), the residual R V (ν) = 0. Then at iteration (ν + 1), we try to correct the solution V (ν) by a correction δV (ν+1) (= V (ν+1) -V (ν) ) calculated by:

∂R (ν) ∂ V δV (ν+1) = -R V (ν) (2.42)
New solution can be then update by

V (ν+1) = V (ν) + δV (ν+1) .
For the resolution of Eq.2.40 with Newton-Raphson method, two approaches exist depending on the primary solution V , correspondingly the enthalpy or temperature. The resolution based on temperature will be presented and used in the following. Thus the system to be resolved is:

∂R ∂ T (ν) ij T (ν+1) j -T (ν) j = -R (ν) i (2.43)
The most important is to compute the tangent matrix ∂R ∂ T ij , which controls the convergence speed. It can be easily deduced from Eq.2.40:

∂R ∂ T (ν) ij = M ij ∂H ∂ T (ν) j no sum on j + A ij ∂H ∂ T (ν) j no sum on j +K1 ij + K2 ij (2.44)
The computation of tangent matrix requires the evaluation of ∂H/∂T at each node j. This can be done by the formulation of Morgan [START_REF] Desmaison | Modélisation Numérique d'un Procédé de Soudage Hybride Arc / Laser En Approche Level Set : Application Au Soudage Multi-Passes de Tôles d'acier de Forte Épaisseur[END_REF]:

∂H ∂ T (ν) j = H (ν) j -H (ν-1) j T (ν) j -T (ν-1) j (2.45)
In the case of the first iteration (1) of Newton-Raphson, the above formulation is adapted to:

∂H ∂ T (1) j = H t j -H t-∆t j T t j -T t-∆t j (2.46)
2.6. SINGLE TRACK SIMULATION ON SUBSTRATE

Single track simulation on substrate 2.6.1 Material properties and configuration

As only the heat transfer is modeled at this stage, the following validation cases are based on a laser irradiation on a 100% dense alumina substrate and consequently no displacement of gas/material interface. The fluid dynamics in the melt pool is neglected, thus no convection is considered. The validation cases begin with the heat source model, which is the most important as it is the driving force of all physical phenomena. The implementation of heat source model will be validated by comparison with analytical results. The computation time is also discussed to study the parallelization efficiency. Regarding the temperature distribution and melt pool geometry, the influence of material and process parameters will be investigated.

The energy conservation will also be verified by imposing adiabatic boundary conditions.

Material properties

Material properties used for thermal modeling are essentially density, thermal conductivity, absorption and heat capacity or specific enthalpy. These properties of gas and alumina are given in Table 2.2. Among them, a constant reflection coefficient R of substrate surface is simply taken to be 0.05 (near to the averaged value measured for powder bed surface, detailed later in Section 5.2.1) for this sensitivity investigation. Although the shielding gas is usually argon, as no data have been found for the surface tension between argon and liquid alumina, the surface tension between air and liquid alumina will be used in the modeling of melt pool dynamics in the next chapter. Consequently, all gas properties are taken as those of air hereafter in order to keep the consistence.

Material

Property Symbol Value Unit Ref.

Gas

Density

ρ g 1.3 kg • m -3 [87] Conductivity λ g 0.024 W • m -1 • K -1 [87] Specific heat capacity C g p 1000 J • kg -1 • K -1 [87] Absorption α g 0 m -1 Alumina Density ρ a 3970 kg • m -3 [88] Conductivity λ a Figure 2.15(a) W • m -1 • K -1 [89] Specific enthalpy h a Figure 2.15(b) J • kg -1 [90] Other Reflection R 0.05 Table 2.2: Material properties
The conductivity of alumina is plotted in Figure 2.15(a) according to the following expression [START_REF] Touloukian | Thermal expansion -nonmetallic solids[END_REF]:

λ a = 5.5 + 34.5 exp -3.3 × 10 -3 T T ∈ [25, 1300] • C (2.47)
and the values are extrapolated for T > 1300 • C. Conductivity values are relatively higher than other ceramics like zirconia (∼ 2 W • m -1 • K -1 ), especially at low temperature. The specific enthalpy of solid and liquid alumina can be found from database [START_REF] Chase | Thermalchemical Tables[END_REF] and modeled as a polynomial function of temperature. As an artificial solidification interval from 2004 to 2104 • C is assumed, the specific enthalpy is averaged between solid and liquid in this interval, leading to a jump of specific enthalpy as shown in

Simulation configuration

All validation cases use the same system configuration as illustrated in Figure 2.16, with a total dimension of 3×0.5×1.1 mm 3 . The substrate is at the bottom with a height of 0.95 mm and the rest is gas. A laser scans the powder along the median plan from X S = 0.2 mm to X E = 2.8 mm. Initial and boundary conditions are:

• T 0 = T ext = 20 • C;
• Top -adiabatic;

• Bottom and lateral faces -heat exchange with

h c = 40 W • m -2 • K -1 .
In order to verify the heat source distribution with non constant absorption, we will assume an absorption field depending on the level set function ψ. An analytical solution of qL can be obtained with which the numerical simulation will be compared. Then the influence of different parameters on the temperature distribution and melt pool shape is investigated. Two groups of study cases are carried out regarding the material properties related to heat source (absorption) and process parameters (interaction radius, laser powder and scanning speed).

The time step of simulation is set to ∆t = 4 µs. Mesh adaptation is used to track the melt pool boundary. Details of adaptation method will be presented in the next chapter. The number of FE element used is 30 000 at the beginning and increased to 50 000 at the end of the simulation.

Heat source distribution

The implementation of heat source model (detailed in Section 2.5.1) is firstly validated. The heat source distribution should be well described even with variable absorption coefficient α. Hence we consider a linear distribution of material absorption expressed by:

α D 1 (ψ) =    10(1 + ψ) if ψ ≤ 0 0 if ψ > 0 (2.48)
Here ψ is in mm and α D 1 in mm -1 . Inserting Eq.2.48 into Eq.2.22, one can easily get the analytical heat source distribution as:

qL (r, ψ) =      (1 -R) • 2P L πr 2 int exp -2r 2 r 2 int • 10(1 + ψ) • exp 10ψ -5ψ 2 if ψ ≤ 0 0 if ψ > 0 (2.

49)

P L =84 W and r int =37.5 µm are used in the simulation. The grid size in (X, Y, Z) directions are (2, 2, 5) µm, respectively. It is finer in X and Y directions in order to better approximate the Gaussian distribution. Figure 2.17(a) shows the heat source distribution. The value of qL along the laser axis (white line) is plotted in Figure 2.17(b) with scatter points. One can see that the simulation matches well the analytical solution. In the transition zone around ψ = 0, there is one point with very small value. This is due to the average of absorption coefficient between gas and material, leading to the global absorption {α} in ψ ∈ -, 0 smaller than that of α D 1 and a smaller qL according to Eq.2.22.

The iso-contours of heat source distribution are also compared in Figure 2.18. The result shows a good coherence between simulation and analytical solution. Note that the precision of simulation depends on the mesh size and also the grid size. Regarding these results, we validate the implementation of heat source model. The computation time is investigated for the case with and without relocalization (search of the containing element of each grid point) as detailed in Table 2.3. The computation time for the localization of grid points is about 89% of the total time in the case with relocalization. However, if the relocalization is not needed, the localization relationship of grid points is taken from memory. This process is almost instantaneous. Note that the relocalization is done only one time between two mesh adaptations.

Influence of material properties

The influence of material properties on the temperature distribution and melt pool shape is firstly investigated. Specifically, the studied material property is the absorption of alumina, in solid and liquid states, as shown in Table 2.4. Firstly, an additional case with the same parameters as case #1 but with adiabatic boundary conditions and without radiation at the material surface is done to verify the energy conservation. The total enthalpy and the total energy received from laser are plotted every 100 time steps (for better visualization) as shown Figure 2.19. The results show that good energy conservation is achieved.

# α s α l r int P L v L mm -1 mm -1 ( µm) (W) mm • s -1
The evolution of temperature field in reference case #1 is shown in Figure 2.20. At the transition stage before the steady state, the melt pool is developed progressively and elongated in the scanning direction. At t = 6 ms, the melt pool shape is almost not changed compared with t = 4 ms, similarly for the iso-contours T = 2004 and 1500 • C. However, the iso-contour T = 1000 • C is furtherly elongated. Due to the scanning direction, the thermal gradient in front of the laser is much higher than that in tail of the melt pool. steady stage at t = 6 ms. The increase of absorption to 10 mm -1 (both for solid and liquid) in case #2 leads to a longer and larger melt pool, similarly for the temperature iso-contours. However, the melt pool becomes shallower. Although the local absorption is higher in deep region in case #2 than in #1, the attenuation during the penetration in material is also higher. This results in a lower qL in deep region in case #2, thus a shallower melt pool. Case #3 with high absorption for solid (same as case#2) and low value for liquid (same as case #1) is somehow like an average between the case #1 and #2. The melt pool or iso-contours on the material surface are more elongated than in case #1 and less elongated than in case #2. However, this is not the case in depth, which shows different tendency for different isocontours. For example, the iso-contours T = 2104 and 2500 • C are deeper than those in case #1, while the one of T = 1000 • C is shallower. The result shows a longer, larger and deeper melt pool than that in case #1. In contrast to case #3, case #4 is with high liquid absorption (same as case #2) and low solid absorption (same as case #1). The melt pool length and width in this case is between the case #1 and #2, while it is shallower than that in case #2 and case #1.

Influence of process parameters

The important influence of material absorption is demonstrated in the previous section. However, it is usually determined by material and can not be modified to get suitable melt pool shape. By contrast, process parameters are fully controlled by user. Consequently, it is important to study the influence of process parameters on the melt pool shape and establish the process window. Here the influence of laser radius, power and scanning speed are investigated as shown in Table 2 Case #5 with increased interaction radius decreases energy concentration. The melt pool shape and temperature iso-contours (Figure 2.24) on material surface are almost the same as those in case #1 (Figure 2.21). The main influence is noted in the depth (Figure 2.25).

.5. # α s α l r int P L v L mm -1 mm -1 ( µm) (W) mm • s -1
As mentioned for Figure 2.10, iso-values of qL becomes shallower with increased r int , hence melt pool and iso-temperatures are always shallower than those in case #1. The linear energy

#5 #6 #7 T ( • C) X Y Figure 2
.24: Top view of temperature distribution at the material surface with temperature iso-contours corresponding to T = 1000, 1500, 2004 (solidus), 2104 (liquidus, red), 2500, 4000 and 6000 • C at t = 13 ms, for cases #5-7.

E l = P L /v L can be defined to study the effect of laser power and scanning speed. It represents the energy deposition per scanning length. The increase of power leads to higher E l , hence larger, wider and deeper melt pool and temperature iso-contours. By contrast, the increase of scanning speed decreases E l , leading to smaller melt pool and iso-contours. These are clearly shown in Figure 2.24 and Figure 2.25. Figure 2.26 compares the melt pool geometries of cases #1 and 5-7. The increase of laser power has the same enlarging effect on the length, width and depth, opposite for the effect of scanning speed.

CHAPTER 2. HEAT TRANSFER

The influence of material properties and process parameters on the melt pool dimension is concluded in Table 2.6, together with maximum temperature attended in each case. Among them, the absorption coefficient is crucial in all 3 dimensions, while the interaction radius has more effect on the depth. The influence of laser power and scanning speed are both very significant. Note that the reflection coefficient has opposite influence compared with the laser power as the net energy input is P L (1 -R). For the maximum temperature, higher T max is obtained with increased absorption, especially that of liquid. For process parameters, higher laser power leads to higher T max , while it is in contrast by increasing interaction radius or scanning speed.

# α s α l r int P L v L Length Width Depth T max mm -1 mm -1 ( µm) (W) mm • s -1 (mm) (mm) (mm) ( • C) 1 5 5
2 ↑ ↑ ↑ ↑ ↓ 7623 3 ↑ ↑ ↑ ↑ 4778 4 ↑ ↑ ↑ ↓ 6515 5 ↑ ↓ ↓ ↓ 3187 6 ↑ ↑ ↑ ↑ 6200 7 ↑ ↓ ↓ ↓ 3087
Table 2.6: Influence of investigated parameters modified with respect to reference case #1.

Conclusion

After a review of modeling scales and methods for powder bed based AM processes. The SLM process is modeled at the track scale with FE and level set method. Powder is considered as a continuous medium and particles are not represented, thus reducing the computational charge. Different heat source models are discussed. The model based on the Beer-Lambert law taking into account the material absorption is employed and implemented in the context of continuous powder and weak absorption of ceramic materials. The application of this model to variable absorption coefficient is validated. Several validation tests are carried out to investigated the influence of material properties and process parameters on the melt pool shape. The significant influence of absorption coefficient is pointed out. In addition, the distribution of absorption can largely affect the heat distribution and thus the melt pool shape as in the cases #3 and #4 (Figure 2.21-2.23). Process parameters like laser radius, power and scanning speed can also largely change the melt pool. Consequently, in the real process, we can adapt process parameters in order to obtain suitable melt pool with desired dimensions.

Chapter 3

Fluid dynamics in melt pool

Résumé

Avec l'hypothèse des milieux continus pour le lit de poudre, la modélisation de la dynamique du bain liquide est réalisée en couplant les résolutions des équations de conservation de la masse et de la quantité de mouvement, et en utilisant une loi de comportement newtonienne compressible pour la poudre et newtonienne incompressible pour le liquide et l'air. L'évolution de la surface du bain liquide est suivie par la level set ψ = 0. La modélisation et l'implémentation des forces motrices est nécessaire pour la dynamique du bain liquide. Dans le modèle, la gravité constante, la tension de surface et la force de Marangoni sont intégrées. Malgré l'importance de l'effet de la vaporisation, à travers la pression de recul, ces dernières ne sont pas encore intégrées dans notre modèle. L'implémentation du modèle par éléments finis est ensuite présentée, notamment la formulation semi-implicite de la tension de surface. L'adaptation du maillage est également détaillée, vue son importance pour le suivi de l'interface mobile.

Le modèle est d'abord validé par un cas simple. Ensuite, l'influence de la tension de surface, de la viscosité et de la vitesse de balayage laser sur la dynamique du bain liquide est étudiée par simulation monopasse du lasage d'une couche de poudre sur un substrat. L'effet Marangoni modifie le flux de convection à l'intérieur du bain liquide, ainsi que sa forme, avec différentes influences dépendant du signe du coefficient de Marangoni ∂γ/∂T . Cet effet redistribue la chaleur et diminue la température du bain liquide. De plus, le signe de ∂γ/∂T a une influence directe sur la morphologie de la surface du cordon. Une autre propriété importante du matériau est la viscosité du liquide, qui a un effet opposé à celui de la tension de surface au regard de la géométrie de l'interface gaz/matière. Plus la viscosité est élevée, moins la rugosité de surface est élevée. Les calculs montrent que la simulation développée est capable de rendre compte de ces divers effets et de détecter l'apparition de l'effet de "balling", où la qualité de la surface du cordon est dégradée lorsque la vitesse atteint certaines limites.

STATE OF THE ART

This chapter focuses on the modeling of fluid dynamics in melt pool. By a brief literature review, the modeling of different physical phenomena contributing to the melt pool dynamics is presented, such as surface tension, Marangoni effect and recoil pressure (induced by evaporation). The modeling method [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramic materials[END_REF] in this work is then detailed, with a compressible Newtonian behavior for the assumed continuous powder and incompressible Newtonian behavior for the liquid. Tracking of the melt pool surface is conducted by the level set method. Driving forces taken into account in this model are gravity, surface tension and Marangoni force. The numerical implementation of these forces with the finite element method is presented, together with mesh adaptation for tracking the evolving gas/material interface. At the end, several cases are presented to demonstrate the capacity of the model and investigate the influence of surface tension, viscosity and scanning speed on melt pool dynamics.

State of the art

Multiple physical phenomena are coupled together in powder based AM processes. The thermal phenomena (including laser/material interaction and heat transfer) melt powder and result in the formation of melt pool. Driving forces acting on the melt pool include gravity, buoyancy, surface tension, Marangoni force and recoil pressure induced by evaporation (see Section 1.3.2), as presented by Markl et al. [START_REF] Markl | Multiscale Modeling of Powder Bed-Based Additive Manufacturing[END_REF]. The resulting dynamics in melt pool is essential to the heat and mass transfer. It changes the melt pool and track shape, the distribution of chemical components and thus the solidification structure. High dynamics may cause the well known balling effect and the spattering of liquid droplets, which are undesirable. Experimental observations and in-situ measurements of melt pool dynamics are difficult on such a small scale. Numerical simulation is expected to bring light on the melt pool dynamics and help to establish process windows.

Surface tension is non negligible in the modeling of melt pool dynamics. As introduced in Section 1.3.2, it can smooth the melt pool surface or lead to balling effect, depending on the scanning velocity. The effect of surface tension was studied by Khairallah and Anderson [START_REF] Khairallah | Mesoscopic simulation model of selective laser melting of stainless steel powder[END_REF] as shown in Figure 3.1. Simulation results suggest that with surface tension, molten particles are bound together and form a melt pool with smoother surface, while without surface tension, molten particles are almost isolated. The binding effect with surface tension can increase the contact and favor heat transfer. King et al. [START_REF] King | Overview of modelling and simulation of metal powder bed fusion process at Lawrence Livermore National Laboratory[END_REF] investigated the surface tension effect combining the influence of scanning speed. In Figure 3.2 with (a) high scanning speed on a overhang region, severe balling occurs with fragmentation of the melt pool. By (b) decreasing the scanning speed while keeping almost the same linear energy, more continuous track is achieved. The track in the second layer in Figure 3.2(c) is continuous as the (b) previous consolidated layer (serving as substrate) helps stabilizing the top melt track.

The Marangoni effect caused by the surface gradient of the surface tension induces tangential convection flow. As the surface tension is usually temperature dependent, the high temperature gradient between the region under laser spot and lateral sides or the rear of the melt pool results in high melt pool flow. This can significantly change the temperature distribution and the melt pool shape. Yuan and Gu [START_REF] Yuan | Molten pool behaviour and its physical mechanism during selective laser melting of TiC/AlSi10Mg nanocomposites: Simulation and experiments[END_REF] developed a 3D finite volume model for TiC/AlSi10Mg considering the Marangoni effect. A wider and shallower melt pool was observed when this effect was taken into account, giving a melt pool shape close to experiment, as shown in Figure 3.3. However, due to the absence of track shape above the substrate, the influence of Marangoni effect on the track shape evolution could not be investigated. Khairallah et al. [START_REF] Khairallah | Laser powder-bed fusion additive manufacturing: Physics of complex melt flow and formation mechanisms of pores, spatter, and denudation zones[END_REF] improved the model in [START_REF] Khairallah | Mesoscopic simulation model of selective laser melting of stainless steel powder[END_REF] by the addition of Marangoni effect and recoil pressure. They reveal the generation of pore defects, spattering and denudation induced by strong melt pool flow. Figure 3.4(a) shows the melt flow with constant surface tension. Melt agglomeration is formed to minimize surface energy driven by surface tension. By taking into account the Marangoni effect by temperature dependent surface tension in Figure 3.4(b), melt flow is created from hot region (red) under laser spot to colder ones at the rear (green). This convection flow takes heat away from the hot region, hence decreases the maximum temperature. Ejection of liquid droplets (or spattering) with low viscosity is also observed in the simulation. The integration of recoil pressure to account for the evaporation phenomenon in Figure 3.4(c) adds an extra force on the melt pool surface and leads to its depression. The melt pool depth is significantly increased. The keyhole effect may be reached for high recoil pressure.

The surface depression can also result in the formation of porosity as shown in Figure 3.5. The recoil pressure (a) pushes the liquid melt away from the laser impact zone and forms a concave melt pool surface, (b) the lateral liquid tends to fill the cavity under the effect of gravity and surface tension and (c) the liquid along the lateral sides slides and connects together, thus entrapping the gas. If the gas can not escape before solidification, porosity is then generated in the final part. The modeling of keyhole and the formation of porosity caused by recoil pressure can be also found in welding. Courtois et al. [START_REF] Courtois | A new approach to compute multi-reflections of laser beam in a keyhole for heat transfer and fluid flow modelling in laser welding[END_REF] presented a detailed model of recoil pressure. In addition, the multi-reflections of laser beam in the keyhole is modeled using Maxwell's equations. With the commercial code Comsol, they are able to simulate the high deformation of melt pool surface and bubbles during the collapse of the keyhole. These bubbles lead to the residual porosity if they can not escape before the solidification.

Although the simulations with powder representation can predict more details and more physical phenomena, they are very time-consuming. As mentioned by Khairallah and Anderson [START_REF] Khairallah | Mesoscopic simulation model of selective laser melting of stainless steel powder[END_REF], their simultions consume on the order of 100 000 cpu hours. On the other hand, simulations without track shape is less precise. In the following, the fluid dynamics will be modeled. The modeling at track scale with powder representation requires less computation resource while it can still give access to the track shape.

Modeling of fluid dynamics in melt pool

Governing equations

The modeling of fluid dynamics is established by coupling the mass and momentum (Navier-Stokes) conservation equations. Due to the change of apparent density during melting from powder to liquid, compressible flow must be taken into account:

         ∇ • u = θ ρ ∂u ∂ t + ∇u • u -∇ • σ = f v (3.1)
with boundary conditions:

     u = u imp on Γ u σ • n Ω = f imp on ∂Ω f (3.2)
where u is the velocity, σ the stress tensor, f v the volumetric force and θ is the volume variation rate related to the shrinkage from powder to dense matter. n Ω is the outward normal direction along the system boundary. u imp and f imp are given velocity and surface force imposed on boundary ∂Ω u and ∂Ω f , respectively. Although we are only interested in fluid dynamics in melt pool, these two equations are applied to the whole system. Here all variables are in final form averaged between the two domains (material and gas) while the brace { * } will be omitted by default for simplicity if not specified. Note that compressibility is taken into account by θ = 0 and the compressible behavior relating σ to the velocity field. These will be detailed in the following.

Mass conservation

The condensation of material from powder state (seen as continuous) to dense matter leads to the change of its apparent density. Mass conservation should be respected as no mass flux is introduced into the system. The mass conservation is first established both in the material (D 1 ) and the gas (D 2 ) domain:

         ∂ ρ D 1 ∂ t + ∇ • ρu D 1 = 0 ∂ ρ D 2 ∂ t + ∇ • ρu D 2 = 0 (3.3)
where the * bracket stands for the volume averaging in a multiphase medium (see Section 2.2.2.2). In our case, the gas domain (D 2 ) is assumed monophasic, so that ρ D 2 = ρ g and ρu D 2 = ρ g u. In the classical velocity-pressure formulation developed for the resolution of the momentum conservation equation, the divergence of velocity, ∇•{u}, is required. In order to deduce it from Eq.3.3, we can firstly show:

∇ • u D i = - 1 ρ D i   ∂ ρ D i ∂ t + ∇ ρ D i • u D i   ∀ i ∈ {1, 2} (3.4)
with the assumption that ρu D i = ρ u D i . In fact, in the case of constant gas density, ∇ • u D 2 = 0. This leads to the global mass conservation equation:

∇ • u = ∇ • {u} = ∇ • (1 -H) u D 1 + H u D 2 = (1 -H)∇ • u D 1 + H ∇ • u D 2 =0 +∇H • u D 2 -u D 1 = - 1 -H ρ D 1   ∂ ρ D 1 ∂ t + ∇ ρ D 1 • u D 1   + δ u D 2 -u D 1 • n (3.5)
where ∇H = dH dψ ∇ψ = δn and n is the unit normal direction outward of material domain (Eq.2.6) and δ is the smooth Dirac function (Eq.2.11).

Momentum conservation

For the fluid mechanics, a simple Newtonian behavior law is used for both gas and material, at any temperature. For the material, it is considered as compressible in the condensation stage from powder to dense matter, while it is incompressible in the powder and liquid state. In the solidification temperature range from 2004 • C to 2104 • C (mushy zone), the Darcy term is not considered, but the viscosity is considered to increase with the solid fraction and it attains a high fixed value for complete solid (more details in Section 3.6.1). On the other hand, gas is always taken as incompressible. Both compressible or incompressible Newtonian behaviors are addressed by the following equation:

σ = s -pI = 2η ˙ - 1 3 tr ˙ I -pI (3.6) and tr ˙ = ∇ • u = θ (3.7)
where s is the deviatoric part of σ, ˙ the strain rate tensor, I the identity tensor, p the pressure and η the dynamic viscosity. Note that if tr ˙ = 0, Eq.3.6 is simplified to the incompressible law. Inserting Eq.3.6 and 3.7 into Eq.3.1, the Navier-Stokes equation with compressible behavior can be obtained:

ρ ∂u ∂ t + ∇u • u -2∇ • η ˙ + ∇p + 2 3 η∇(∇ • u) = f v (3.8)
With the assumption of constant η, it can be written as:

ρ ∂u ∂ t + ∇u • u -2η∇ • ˙ + ∇p + 2 3 η∇(∇ • u) = f v (3.9)

Integrated driving forces

Several forces contribute to the melt pool dynamics. Among them, the importance of surface tension (including Marangoni effect), buoyancy force and recoil pressure is well recognized. Among them, the buoyancy is relatively less significant than others and the recoil pressure is less critical for ceramics due to their high boiling temperature (e.g. 3240 K for Al 2 O 3 ). Hence, only gravity, surface tension and Marangoni force have been implemented in this work.

Surface tension

The phenomenon of surface tension is caused by the cohesive forces between liquid molecules [START_REF] Bush | MIT Lecture Notes on Surface Tension[END_REF]. For each molecule in the bulk of the liquid, all cohesive forces from neighboring liquid molecules are in local equilibrium. However, at the liquid/gas interface, these forces are partly contributed by gas molecules. As the intermolecular distance in gas is higher than that in

MODELING OF FLUID DYNAMICS IN MELT POOL

liquid, the cohesive forces induced by gas molecules are smaller than those induced by liquid molecules. Consequently, internal pressure is created and the interface has tendency to bend towards the gas. This is described by the Young-Laplace equation [START_REF] De Gennes | Capillarity and Wetting Phenomena -Drops, Bubbles, Pearls, Waves[END_REF]:

p g -p l = -γ∇ • n l = γκ t (3.10)
where p l , p g are pressures in liquid and gas, respectively. γ is the surface tension coefficient and κ t = -∇ • n l the total curvature. n l is the unit normal direction pointing out of liquid and it can be easily computed by an additional level set function ψ l associated to the melt pool boundary:

n l = ∇ψ l ∇ψ l = 1 (3.11)
The resulting surface tension force acting on the liquid surface is:

f s = γκ t n l (3.12)
This is a surface force acting in the normal direction toward the liquid at the liquid/gas interface. The tangential effect of surface tension (or Marangoni effect) will be presented in the following.

Marangoni effect

The importance of Marangoni effect is well known in welding [2,[START_REF] Desmaison | Modélisation Numérique d'un Procédé de Soudage Hybride Arc / Laser En Approche Level Set : Application Au Soudage Multi-Passes de Tôles d'acier de Forte Épaisseur[END_REF][START_REF] Mills | Marangoni effects in welding[END_REF] and SLM [START_REF] Wang | Mechanisms and characteristics of spatter generation in SLM processing and its effect on the properties[END_REF][START_REF] Olakanmi | Selective laser sintering/melting (SLS/SLM) of pure Al, Al-Mg, and Al-Si powders: Effect of processing conditions and powder properties[END_REF]. This effect is due to the shear stress at the liquid/gas interface induced by the gradient of surface tension coefficient along this interface. Assuming that γ essentially depends on temperature, we have:

f M = ∇ s γ = ∂γ ∂ T ∇ s T (3.13)
where ∇ s denotes the surface gradient operator. For instance, ∇ s T is the tangential component of the gradient at the interface:

∇ s T = ∇T -(∇T • n l )n l (3.14)
According to Eq.3.13, with a given surface temperature gradient, the Marangoni force depends only on the Marangoni coefficient ∂γ/∂T , which is a property of the material in liquid state. The sign of ∂γ/∂T plays an important role as the direction of the convection flow depends on it as shown in Figure 3.6. Consequently, temperature distribution and melt pool shape [START_REF] Mills | Marangoni effects in welding[END_REF] will be modified.

The importance of Marangoni effect can be evaluated by the non-dimensional Marangoni number M a [START_REF] Maroto | Introductory analysis of Bénard-Marangoni convection[END_REF], which is the ratio between thermal surface tension forces and viscous forces, expressed by:

M a = - ∂γ ∂ T L∆T ηα T (3.15)
where L is the characteristic length and α T = λ/ρC p the thermal diffusivity. Note that both the thermal conductivity and viscosity can stabilize the Marangoni flow. The former is preemptive as it decreases the temperature gradient, thus the Marangoni force, while the latter slows down the resulting flow. As the viscosity of ceramics is much higher than metals (see Section 3.6.1), the resulting Marangoni flow is certainly less important in liquid ceramics than in liquid metals.

Numerical resolution

Semi-implicit formulation of surface tension

According to Eq.3.12, the surface tension depends on the total curvature κ t and the normal direction n l . However, the interface geometry is not known a priori. In fact, it should be determined by the previous one and the current velocity to be solved. A simple explicit formulation consists in computing κ t and n l based on the interface at previous time increment. However, this formulation is limited by the time step [START_REF] Hysing | A new implicit surface tension implementation for interfacial flows Modeling of Surface Tension Effects[END_REF]:

∆t < l 3 2 ρ 2πγ (3.16)
where l is the mesh size. A simple calculation for alumina and a typical mesh size for simulation of SLM with ρ = 3970 kg • m -3 , γ = 0.64 N • m -1 and l = 1 µm gives a time step smaller than 31 ns, leading to a non-sustainable computational charge. This limitation can be highly reduced by a semi-implicit formulation of surface tension.

The total curvature and normal direction at the liquid/gas interface can be related to the interface position by the following equation [2,[START_REF] Hysing | A new implicit surface tension implementation for interfacial flows Modeling of Surface Tension Effects[END_REF][START_REF] Khalloufi | High fidelity anisotropic adaptive variational multiscale method for multiphase flows with surface tension[END_REF]:

κ t n l = ∆ s x = ∇ s • (∇ s x) (3.17)
where x ∈ Γ denotes only the points at liquid/gas interface Γ. ∆ s and ∇ s (see Eq.3.14 for the application to a scalar) are the surface Laplace and Nabla operators, respectively. When 66

NUMERICAL RESOLUTION

∇ s is applied to a vector x, it can be related to ∇x by:

∇ s x = ∇x -(∇x • n l ) ⊗ n l = ∂ j x i -n j n k ∂ k x i e i ⊗ e j (3.18)
where e i and e j are basis vectors and the index l of n l is omitted in order to avoid ambiguity. Using x = x -+ u∆t, Eq.3.17 at current time step can be expressed by:

κ t n l = ∆ s x = ∆ s x -+ ∆ s u∆t = (κ t n l ) -+ ∆ s u∆t (3.19)
Here the current time step t is omitted for simplicity and the sign minus represents the previous time step. Compared with the explicit formulation, the second term at RHS related to the surface diffusion of velocity is added in this semi-implicit formulation. Although the velocity is taken at current time, the normal direction in ∇ s operator is still taken at the previous time step. its dot product with a vector test function w [START_REF] Hysing | A new implicit surface tension implementation for interfacial flows Modeling of Surface Tension Effects[END_REF]:

Γ γ(κ t n l ) -• w dΓ + ∆t Γ γ∆ s u • w dΓ = Γ γ(κ t n l ) -• w dΓ + ∆t ζ γ(∇ s u • w) • n ζ dζ -∆t Γ γ∇ s u : ∇ s w dΓ (3.20)
where ζ = ∂Γ is the contour line (or boundary of wetted region) of the liquid/gas interface Γ.

Here the integration by part is applied to get the last two terms at RHS. The second term at RHS will vanish if the interface Γ is closed or by taking a test function vanishing at ζ. The contraction part in the last term at RHS can be related to normal Nabla operator [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramic materials[END_REF]:

∇ s u : ∇ s w = ∂ j u i -n j n k ∂ k u i e i ⊗ e j : ∂ j w i -n j n l ∂ l w i e i ⊗ e j = ∂ j u i ∂ j w i -∂ j u i • n j n l ∂ l w i -n j n k ∂ k u i • ∂ j w i + n j n k ∂ k u i • n j n l ∂ l w i n j n j =1 = ∂ j u i ∂ j w i -n j n k ∂ k u i • ∂ j w i = ∇u : ∇w -(∇u • n l )(∇w • n l ) (3.21)
where the index l of n l is also omitted.

Finally, in the level set framework, the surface tension in Eq.3.20 can be transformed into volume force [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] by multiplying the Dirac function δ l , which is related to ψ l and only applied at liquid/gas (part of ψ) or liquid/powder interface (see blue and white curve in Figure 2.7). With Eq.3.21, it becomes:

Ω δ l γ (κ t n l ) -• w -∆t ∇u : ∇w -∇u • n - l • ∇w • n - l dΩ (3.22)
It should be noted that the normal direction n l is taken at previous time step and correspondingly the total curvature κ t .

FE implementation

The resolution of Eq.3.1 is carried out in the framework of a mixed velocity-pressure FE formulation. In the FE context, the choice of functional spaces for velocity and pressure is essential to the resolution stability. The classical Galerkin method with linear P 1 formulation for both velocity and pressure is not stable as it does not satisfy the inf-sup stability condition [START_REF] Babuška | Error-bounds for Finite Element Method[END_REF][START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers[END_REF]. Consequently, oscillation of the solution is encountered. This difficulty can be overcome by two methods: i) the stable mixed finite elements or ii) the Variational MultiScale (VMS) method. The former uses P 1 + /P 1 formulation with enriched degree of freedom for velocity by adding a bubble function [START_REF] Arnold | A stable finite element for the stokes equations[END_REF]. Details of this method is not discussed here but will be introduced in the next chapter when describing the solid mechanic solver. The latter (VMS) was originally proposed by Hughes [START_REF] Hughes | Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods[END_REF] and will be briefly presented in the following. The solution is considered at two scales: the coarse and fine scales. Both the velocity and pressure are decomposed into these two scales. The 3.3. NUMERICAL RESOLUTION resolution of Eq.3.1 can be also split into two scales. The solution at fine scale is then used in the resolution at coarse scale, giving the final solution. More details of treatment can be found in the work of Hachem [START_REF] Hachem | Stabilized Finite Element Method for Heat Transfer and Turbulent Fows inside Industrial Furnaces[END_REF] and Khalloufi et al. [START_REF] Khalloufi | High fidelity anisotropic adaptive variational multiscale method for multiphase flows with surface tension[END_REF] at Cemef. Here we just present the final equations to be solved with semi-implicit surface tension and compressible behavior.

Firstly we can define the following functional spaces:

V = u | u ∈ H 1 (Ω) d , u = u imp on ∂Ω u (3.23) V 0 = u | u ∈ H 1 (Ω) d , u = 0 on ∂Ω u (3.24) P = p | p ∈ L 2 (Ω) (3.25)
The resolution consists in finding u ∈ V and p ∈ P, such that ∀w ∈ V 0 and q ∈ P:

                               ρ u ∆t , w supg + ρ∇u • u, w supg + 2 η ˙ (u) : ˙ (w) + ∇p, w supg + (τ C ∇ • u, ∇ • w)-2 3 (η∇ • u, ∇ • w) +∆t(δ l γ∇u : ∇w) -∆t δ l γ∇u • n - l • ∇w • n - l = f t , w supg + τ C θ, ∇ • w (∇ • u, q) + ρ u ∆t , τ K ∇q + (ρ∇u • u, τ K ∇q) + (∇p, τ K ∇q) = θ, q + (f t , τ K ∇q) (3.26) with w supg = w + τ K ρ∇w • u (3. 27
)
f t = ρg + δ l (κ t n l ) -+ ∂γ ∂ T ∇T -(∇T • n l )n l + ρ u - ∆t (3.28)
Here the notations ( * , * ) = Ω ( * • * ) dΩ and ( * : * ) = Ω ( * : * ) dΩ are used. Terms in blue in Eq.3.26 are the contribution of semi-implicit surface tension (the explicit contribution is hidden in f t ) while the ones in red are related to compressibility. The stabilization coefficients τ k , τ C are taken as the same in [START_REF] Khalloufi | High fidelity anisotropic adaptive variational multiscale method for multiphase flows with surface tension[END_REF][START_REF] Coupez | Solution of high-Reynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing -ScienceDirect[END_REF]:

τ K =   2ρ ∆t 2 + 2ρ u h K 2 + 4η h 2 K 2   -1/2 (3.29) τ C =   η ρ 2 + c 2 u c 1 h K 2   1/2 (3.30)
where c 1 and c 2 are constants independent to characteristic length h K . The calculation of h K is proposed Tezduyar and Sathe [START_REF] Tezduyar | Stabilization parameters in SUPG and PSPG formulations[END_REF]:

h K = 2 u i=D i=1 |u • ∇N i | (3.31)
where N i is the interpolation function of node i in the local element.

Remeshing strategy

Mesh quality is essential in our FE method. A suitable mesh should be fine enough to ensure a precise resolution. On the other hand, it should be as coarse as possible in order to reduce the computational charge. In addition, the use of the LS method requires a refinement of the evolving gas/material interface. For the mesh type, anisotropic mesh is preferred as it allows refining only in the direction with large gradient variation. Such an anisotropic mesh adaptation is available in the library Cimlib and will be detailed in the following. An anisotropic mesh is controlled by a metric tensor field M defined at each node i in the basis (e 1 , e 2 , e 3 ).

M i(e 1 ,e 2 ,e 3 ) =       h 11 h 12 h 13 h 12 h 22 h 23 h 13 h 23 h 33       (3.32)
This is a symmetric tensor and it can also be expressed in its principal orthogonal basis (m 1 , m 2 , m 3 ):

M i(m 1 ,m 2 ,m 3 ) =        1 h 2 m 1 0 0 0 1 h 2 m 2 0 0 0 1 h 2 m 3        = 3 k=1 1 h 2 n k m k ⊗ m k (3.33)
This metric is defined in the way that the distance of the edge x ij connecting the node i and j measured in this metric should be equal to 1:

x ij M i = x T ij M i x ij = 1 (3.34)
The geometrical interpretation of Eq.3.34 in 3D case is that all nodes connected to node i is located at an ellipsoid surface with semi-principal axis h m 1 , h m 2 and h m 3 along in direction m 1 , m 2 and m 3 , respectively (see Figure 3.8 for the demonstration of 2D case).

i e 1 e 2 m 1 m 2 x ij j h m 1 h m 2

Metric construction

The metric construction method used hereafter is based on edge length distribution tensor and edge error estimation, presented and implemented in the library Cimlib by Coupez [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF].

The idea is to deduce a metric tensor at each node according to the edge length distribution tensor. The latter is calculated by the desired lengths of edges connecting to this node. The desired edge length can be obtained by analyzing the interpolation error along each edge and then minimizing the total interpolation error. Readers are encouraged to review the article of Coupez [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF] for more theoretical demonstrations.

If all the edges connected to the node i satisfy the unit length condition in the metric space M i (Eq.3.34), they should be located in en ellipse (2D) or ellipsoid (3D). However, the desired edge lengths do not necessarily meet this condition. Hence, an optimal metric at node i can be obtained by minimizing the following function:

M i = arg min M j∈Γ(i) X T ij MX ij -1 2 (3.35)
where X ij is the desired edge vector connecting node i and j. The continuous metric field defined at the mesh node i is then:

M i = 1 d X -1 i and X i = 1 Γ(i) j∈Γ(i) s 2 ij x ij ⊗ x ij (3.36)
where x ij is the actual edge vector and X i is the distribution tensor of desired edge length. Γ(i) the cardinality of the set Γ(i) containing all nodes connecting to node i and s ij = X ij / x ij is the stretching factor. It is defined by the total error e and the local edge error e ij :

s ij = min    e e ij 1 p+2 , x ij h min    (3.37)
where p is a coefficient which is taken equal to 1.5 and h min is the minimal mesh size set by user. The total interpolation error e is calculated by:

e =        NbN i=1 j∈Γ(i) e p p+2 ij D(D -1)NbE        p+2 p (3.38)
Note that the sum is firstly done with all edges connected to node i and then with all nodes. Consequently, each edge is counted two times and the total number is close to D(D -1)NbE.

The only variable to be determined is the edge error e ij and it is estimated by:

e ij = max |G ij • x ij |, err x ij 2 (3.39)
where err is a coefficient set by user to control the global mesh size and its effect will be shown at the end of this section. One can see that this coefficient is equivalent to the norm of the Hessian matrix (second derivative). G ij = G j -G i is the gradient difference with G i the gradient defined at the node i, which will be detailed below.

(∇f

) E k (∇f ) E l x ij i j
Figure 3.9: The gradients of variable in two neighboring elements E k and E l are constant in each elements, but their projections on the common edge x ij are the same.

Considering a P 1 scalar field f , then its gradient ∇f is P 0 (Fig. 3.9), which is constant in the same element and discontinuous from element E k to element E l . However, the projection of two neighboring elements on the common edge vector x ij is continuous because:

∇f • x ij = f j -f i = f ij (3.40)
Then the objective is to find a gradient field at node i which minimizes the total error of edges connected to this node. It corresponds to minimize the following function:

G i = arg min G    j∈Γ(i) (G -∇f ) • x ij 2    = arg min G    j∈Γ(i) (G • x ij -f ij 2    (3.41)
Taking the derivative (with respect to G) of Eq.3.41 equal to 0, we can get:

G i =    j∈Γ(i) x ij ⊗ x ij    -1    j∈Γ(i) x ij f ij    (3.42)
The tracked field f is not limited to one variable. Several fields of interest can be simultaneously used to deduce the metric tensor. In this case, each field has its corresponding gradient field and the edge error G ij • x ij will be summed in Eq.3.39. However, due to the different spatial evolution of different fields, e.g. Heaviside (to track the interface) and liquid fraction (to track the melt pool boundary), their second derivatives may have totally different magnitudes. Hence, treatment like mapping the liquid fraction to a Heaviside function can be employed in order to better represent the melt pool boundary.

In this method, one important parameter is err , which should be carefully chosen according to the global mesh size and that in the zone of interest. Figure 3.10 shows the effect of this parameter in a 2D case with an inclusion in unit square. A Heaviside function is defined to track the interface and evolves continuously in a thickness of 0.01 around the interface. An initial isotropic mesh with mesh size 0.1 is used. The required element number and minimum mesh size are set to be 10000 and 5 × 10 -4 , respectively. The results in Figure 3.10 show that the refinement based on tracked variable (here the Heaviside function) becomes more obvious by decreasing err . However, the zone out of interest has very coarse elements. By contrast, high err results in smaller global mesh size, while the representation of interface is degenerated. Consequently, suitable mesh adaptation requires preliminary study of this parameter. In the case with requirement of fine mesh in the zone of interest while avoiding a too coarse global mesh, the method of mesh intersection can be used, which will be introduced in the following section.

Intersection of metrics

The metric generated by error estimation and length distribution tensor can well refine the region of interest while keeping almost the same total number of elements. However, at the same time, it leads to the coarsening of region with small gradient variation of tracked variable. This degeneration is difficult to control and may sometimes result in a loss of precision. Although one can add tracked variables with high gradient variation in this region to avoid this problem, the combination of several tracked variables is difficult to adjust the weight of each variable as they have different order of magnitude of gradient variation. Thus, It is sometimes useful to employ additional metrics to control the mesh size. For example, a constant metric may be used to control the global mesh size in order to avoid too coarse elements.

The mesh generator accepts only one metric. Therefore, the use of several metrics requires their intersection. The idea is to find the metric with corresponding ellipse (or ellipsoid in 3D) which is inscribed in all ellipses (or ellipsoids in 3D) represented by these metrics. It has been shown in Eq.3.33 and Figure 3.8 that the ellipse of a metric can be defined by its eigenvectors (direction of semi-principal axis) and corresponding eigenvalues (length of semiprincipal axis). Considering a simple case with an isotropic metric M 1 and an anisotropic one M 2 , the intersection of these two metrics can be simply carried out by taking the minimum eigenvalues of two metrics as they have the same eigenvectors (any direction can be taken as eigenvectors for the isotropic one). However, for a general case with two anisotropic metrics, this does not work as they have different eigenvectors.

M 1 ∩ M 2 (M 1 ∩ M 2 ) ∩ M 3 M 1 M 2 M 3 M 2 ∩ M 3 M 1 ∩ (M 2 ∩ M 3 ) M 1 M 2 M 3
Figure 3.11: The intersection of several metrics is not associative (after [START_REF] Lo | Finite Element Mesh Generation[END_REF]).

One intersection method is called "simultaneous diagonalization" [START_REF] Lo | Finite Element Mesh Generation[END_REF]. One can firstly find a common basis (e 1 , e 2 , e 3 ) (not necessarily orthogonal) in which both M 1 and M 2 are diagonal. This basis can be obtained by the diagonalization of M -1 1 M 2 . Then the two metrics can be simultaneously diagonalized by:

P T M 1 P =       λ 1 0 0 0 λ 2 0 0 0 λ 3       and P T M 2 P =       µ 1 0 0 0 µ 2 0 0 0 µ 3       (3.
M 1 ∩ M 2 = P T -1       min (λ 1 , µ 1 ) 0 0 0 min (λ 1 , µ 1 ) 0 0 0 min (λ 3 , µ 3 )       P -1 (3.44)
It should be mentioned that the intersection process is not associative [START_REF] Lo | Finite Element Mesh Generation[END_REF] as illustrated in Figure 3.11. A precise representation of mobile interface is essential. This requires both a robust resolution of transport equation and adequate reinitialization by geometrical method mentioned in Section 2.2.1. Although the reinitialization [START_REF] Shakoor | An efficient and parallel level set reinitialization method -Application to micromechanics and microstructural evolutions[END_REF] method is easy and stable, it is intrinsically non conservative, as the interface can be slightly changed before and after reinitialization. This problem can be explained by the schematics in Figure 3.12 with unit equilateral triangles. Suppose that the transport velocity is zero at time t-∆t and the level set field ψ is given by values reported at FE nodes in (a). By detecting the element edge with sign change of ψ, one can find the intersection point of ψ = 0 with this edge according to the distance |ψ|. Then the interface ψ = 0 (red line in (b)) can be deduced. Based on this interface, the distance of each node to this interface can be geometrically calculated, thus leading to the reinitiallized ψ in (c). At time t, the transport velocity is still zero, while the interface deduced at this time will be slightly different from the previous increment as ψ(A) is changed.

Blockage of elements

1 2 0 0 -1 2 -1 2 (a) A B C D E ψ t -∆ t = 0 = √ 3 4 = 1 2 (b) √ 3 4 0 0 -1 2 -1 2 (c) A B C D E ψ t = 0 ψ t -∆ t = 0 (d)
The stable solution is achieved when the intersection point attends point A, where

ψ(A) = ψ(B) = ψ(C) = 0 and ψ(E) = ψ(D) = - √ 3/2.
This means that the precision is limited to the element size. However, if we consider the mesh adaptation, this 'anchored' intersection point can pass the element, leading to non-physical propagation of the interface. This will degenerate the solidified track surface and should be avoided. Hence, the blockage of elements is proposed when the track surface is solidified. This treatment prevents exceeding transport over more than the mesh size. Although the interface can still change in the scale of element size, good precision is promised if the interface is well refined. The difference of interface evolution with and without blockage will be illustrated in the next section.

Verification of energy and mass conservation

Duo to interpolation error during mesh adaptation, it is essential to verify energy and mass conservation, especially with complex interface evolution. A simple 1D case with fixed homogeneous mesh is thus designed to verify energy and mass conservation with density variation induced by phase transformation. The whole domain is a cuboid with dimensions 0.5×0.5×4 mm 3 (Figure 3.13(a)). Material is below the interface ψ = 0 located at z = 3 mm, and gas is above the interface. At the bottom, there is a homogeneous volume heat source (magenta) in the zone with dimension 0.5 × 0.5 × 1 mm 3 and a power P L = 12.5 W. In the material domain, the initial phase 1 will be transformed to phase 2 with higher density in temperature range [1100 1200] • C (see material properties in Table 3.1). The volume fraction of phase 1 evolves linearly from 1 to 0 in this temperature range. The time step is ∆t=0.1 ms. Initial and boundary conditions are set to be:

• Uniform initial temperature T 0 = 1080 • C, adiabatic conditions on all surfaces.

• Top -free; bottom -non slipping condition with u = 0; lateral faces -slipping condition with u n =0.

Gas Phase 1 Phase 2 In the mass conservation (Eq.3.5), both the velocities of gas and material are engaged in the second term at RHS. However, the resolution coupling with NS equation gives only one velocity field. For the following tests, the second term at RHS of Eq.3.5 will be firstly neglected. Influence of the mesh size l and LS half thickness on the mass conservation are The energy of the whole system can be calculated by:

ρ (kg•m -3 ) 1.3 5680 6816 C p (J•kg -1 •K -1 ) 600 600 600
κ (W•m -1 •K -1 ) 250 250 250
η (Pa•s) 10 -4 10 -3 10 -3
E = Ω {ρh}dV (3.45)
It is compared with the theoretical energy evolution calculated by E th = E 0 + P L t, where E 0 is the initial energy of the system. The evolution of E and E th in Figure 3.14(a) shows that energy conservation is well respected. Note m 0 as the initial mass of material domain and m the current mass, the mass variation ratio is proposed to be evaluated by ∆m/m 0 , where ∆m = m -m 0 . Evolution of ∆m/m 0 is plotted in Figure 3.14(b). There are 3 regimes: 0 ∼ 0.018 s, 0.018 ∼ 0.04 s and > 0.04 s. It can be verified that the first corresponds to phase transformation below the gas/material interface. The second corresponds to the period when phase transformation happens at the interface. The phase transformation finished in the third regime. After the transformation reaches the transition zone, there is always a tendency of mass increase. Mesh size seems to be not related to this problem while this deviation decreases with , which is reasonable to tend to 0 when = 0.

Mass increase may be related to the absence of the second term at the RHS of Eq.3.5. The velocity plotted in Figure 3.13(c) indicates a linear evolution of velocity in ψ ∈ [-, 0]

∇ • u D 1 = ∂ u z D 1 ∂ z = - 1 ρ D 1   ∂ ρ D 1 ∂ t + ∇(ρ) D 1 • u D 1   ≈ - 1 ρ D 1 ∂ ρ D 1 ∂ t (3.46)
Here the approximation in the expression with density gradient ∇ ρ D 1 is due to the negligible spatial variation of ρ D 1 . Then the constant velocity of gas in the transition zone [-, ], u 0 D 2 , is considered as equal to the material velocity, u z D 1 , at the gas/material interface for the continuity reasons. Thus the two averaged velocities can be expressed by:

       u z D 1 = u 0 D 2 -ψ ρ D 1 ∂ ρ D 1 ∂ t ψ ∈ [-, 0] u z D 1 = u 0 D 2 ψ ∈ [0, ] (3.47)
The second term at the RHS of Eq.3.5 is then estimated by averaging the velocity difference in the transition zone ψ ∈ [-, ]:

u D 2 -u D 1 • n ≈ 1 2 - u z D 2 -u z D 1 dψ = 1 4 ρ D 1 ∂ ρ D 1 ∂ t (3.48)
Finally, with δ function, this term is equal to:

δ u D 2 -u D 1 • n = 1 4 ρ D 1 ∂ ρ D 1 ∂ t cos 2 πψ 2 (3.49)
Test #6 with the same parameters as #2 but with the proposed term shows an evident contribution in the mass conservation. However, it shows a slight decrease of mass in the regime 1. Although this proposition is based on a simple 1D case, it will be applied to the cases hereafter.

SINGLE TRACK SIMULATION WITH POWDER ON SUBSTRATE

Single track simulation with powder on substrate

In this chapter, simulation of SLM is demonstrated with fluid mechanics. Due to the importance of melt pool dynamics, the influence of several parameters like liquid viscosity, surface tension (temperature dependent) and scanning speed will be investigated. In particular, the liquid velocity field, temperature iso-contours in the melt pool and the final track shape will be compared when using different mateiral and process parameters.

Material properties and configuration

Alumina powder with 50% porosity will be considered in the following. Besides the properties of alumina presented in Table 2.2, special model should be used for the thermal conductivity of the powder bed. A review on the effective thermal conductivity in packed bed of monosized spherical particles is given by W. van Antwerpen et al. [START_REF] Van Antwerpen | A review of correlations to model the packing structure and effective thermal conductivity in packed beds of mono-sized spherical particles[END_REF]. A frequently used model is proposed by Zehner and Schlunder [START_REF] Zehner | Thermal conductivity of granular materials at moderate temperature (in German)[END_REF]. According to this model, the effective thermal conductivity of powder, λ Z 1 , can be related to those of gas, λ g , and alumina, λ a , by:

λ Z 1 λ g = 1 -1 -g g Z 1 + 2 1 -g g Z 1 1 -ξB   (1 -ξ)B (1 -ξB) 2 ln 1 ξB - B + 1 2 - B -1 1 -ξB   (3.50)
where ξ = λ g /λ a , B = 1.25

1-g g Z 1 g g Z 1
10/9

and g g Z 1 is the powder porosity. With the thermal conductivity of gas (0.024 W • m -1 • K -1 ) and alumina (Figure 2.15(a)), the thermal conductivity of powder with 50% porosity can be obtained as in Figure 3.15(a). The material conductivity λ D 1 is averaged between powder, λ Z 1 , and dense matter, λ Z 2 (same as λ a in Figure 2. [START_REF] Wohlers | History of additive manufacturing[END_REF] The viscosity of liquid alumina is crucial to melt pool dynamics. Langstaff et al. [START_REF] Langstaff | Aerodynamic levitator furnace for measuring thermophysical properties of refractory liquids[END_REF] measured the value by aerodynamic levitation in a wide temperature range from 2173 K (supercooled) to 3240 K (boiling temperature). By fitting their results with an Arrhenius law, we can get the evolution of viscosity depending on temperature:

η l (T ) = 2.15 × 10 -3 exp 57685 R g T [ Pa • s] (3.51)
where R g = 8.314 J • mol -1 • K -1 is the gas constant. Note that the viscosity of liquid alumina at melting temperature T m (=2054 • C) is about 42.2 mPa • s, which is one order of magnitude higher than liquid metals (e.g. 5.5 mPa • s for liquid Fe [START_REF] Battezzati | The viscosity of liquid metals and alloys[END_REF]). Beside the viscosity of liquid, the values of powder and solid are also necessary in such a global resolution of Eq.3.1. As we are not interested in the velocity field of solid which should be considered as a fixed phase, a high solid viscosity of 1000 Pa • s is imposed. In the artificial solidification range between T s = 2004 • C and T l = 2104 • C, a logarithmetic evolution is assumed from 1000 Pa • s to 39.8 mPa • s (red curve in Figure 3.15(b)). For the powder bed, it is difficult to find an adequate model for its viscosity and it is simply assumed to be 10 Pa • s (∼ the viscosity of honey). In the condensation temperature from 1804 to 2104 • C, the viscosity evolution is also assumed to be logarithmetic (black curve in Figure 3.15(b)). The material viscosity η D 1 is then taken between η Z 1 and η Z 2 by a special law:

η D 1 =    η Z 1 if g Z 1 D 1 > 0 (powder exists) η Z 2 if g Z 1 D 1 = 0 (powder totally melted) (3.52)
This means that once the powder is totally melted, η D 1 follows the value of dense matter (Z 2 ) in Figure 3.15(b) and it is revertible. Otherwise, it takes the value of powder (Z 1 ). This choice is taken rather than an arithmetic or logarithmetic mixing law between η Z 1 and η Z 2 in order to avoid the increase of η D 1 during condensation, which causes convergence problem in the numerical resolution.

For the surface tension, we assume that it is linearly dependent on temperature as follows:

γ = γ m + ∂γ ∂ T (T -T m ) (3.53)
where γ m is the surface tension coefficient at melting temperature and the temperature slope ∂γ/∂T is constant. The value of ∂γ/∂T is measured by different methods in literature, varying between -4.8 × 10 -4 and -6 × 10 -5 N • m -1 • K -1 [START_REF] Paradis | Surface Tension and Viscosity Measurements of Liquid and Undercooled Alumina by Containerless Techniques[END_REF]. Among them, the results of Paradis and Ishikawa [START_REF] Paradis | Surface Tension and Viscosity Measurements of Liquid and Undercooled Alumina by Containerless Techniques[END_REF] with γ m = 0.64 N • m -1 and ∂γ/∂T = -8.2 × 10 -5 N • m -1 • K -1 will be used as nominal values hereafter. Due to the uncertainty of ∂γ/∂T , its influence will be investigated.

The simulated system is similar to the Figure 2.16 but with a powder layer of 50 µm deposited on a alumina substrate. Same to the heat transfer simulation in the previous section, P L = 84 W, r int = 37.5 µm, α s = α l = 5 mm -1 and R = 0.05 are set for all cases in the following simulations. In addition, a constant reflection coefficient R = 0.05 is taken, which is near to the experimentally measured average value (detailed later in Section 5.2.1).

SINGLE TRACK SIMULATION WITH POWDER ON SUBSTRATE

Boundary conditions are the same as previous 1D verification tests. Several cases are carried out in order to study the influence of temperature dependent surface tension, liquid viscosity and scanning speed on the melt pool dynamics. Values of these parameters used in each case are detailed in Table 3.3. Case #1 with nominal values of ∂γ/∂T , η l and v L = 200 mm • s -1 will be considered as reference. Then the Marangoni effect can be studied by comparing cases #1, #2 and #3, where case #2 uses the reported extreme value in literature and case #3 aims to study the sign of ∂γ/∂T . Cases #4 and #5 intend to investigate the influence of liquid viscosity and scanning speed, respectively. The mesh adaptation used in these cases is piloted by several metrics. Detailed strategy will be presented when discuss the mesh result. The number of FE elements evolves from ∼ 1 150 000 to ∼ 1 620 000 from the beginning to the end of simulation. The time step is set to be 2 µs for all cases. In the following, the numerical difficulty related to heat convection is firstly discussed. The reference case #1 is presented in details to investigate the temperature distribution, the mesh adaptation and the melt pool dynamics. The conservation of mass and energy is also verified with the reference case. Then the influence of surface tension, viscosity and scanning speed on melt pool dynamics are presented, characterized by velocity field, temperature distribution, melt pool shape and track surface morphology.

# ∂γ/∂T η l (T ≥ 2104 • C) v L (×10 -5 N • m -1 • K -1 ) ( Pa • s) ( mm • s -1 )

Modification of convection velocity in heat transfer equation

The velocity field obtained by resolution of Eq.3.1 includes not only the velocity in the melt pool but also the one induced by the compressible effect, which transports the gas/material interface. If this velocity is used directly in heat transfer, a non-physical phenomena occur as shown in Figure 3.16(a). A temperature drop is observed around the gas/material interface. It causes problems like non-physical solidification of liquid in the region below the interface. The high viscosity of this region furtherly results in numerical difficulties in the resolution of melt pool dynamics.

This phenomenon is due to the convection velocity from gas to material domain, which transports the low enthalpy of gas to material. As the temperature should be in coherence with enthalpy, very low temperature of material should be attained to meet this low enthalpy.

In order to overcome this problem, based on the treatment of Desmaison [START_REF] Desmaison | Modélisation Numérique d'un Procédé de Soudage Hybride Arc / Laser En Approche Level Set : Application Au Soudage Multi-Passes de Tôles d'acier de Forte Épaisseur[END_REF], it is proposed to modify the convection velocity u th in heat transfer equation (Eq.2.23) by:

u th =      (1 -H) u -u • n ψ n ψ if |ψ| ≤ (1 -H)u if |ψ| > (3.54)
This means that only the tangential component of the velocity is kept in the transition zone ψ ∈ [-, ] and the convection in gas is not considered by multiplying by (1 -H). Applying this modified convection velocity in Eq.2.23 results in Figure 3.16(b). The non-physical low temperature can be avoided. This modification does not correspond to the physics of material but can help to avoid the problem of temperature drop. This will be used in the following simulations.

--

SINGLE TRACK SIMULATION WITH POWDER ON SUBSTRATE

Results of reference case

Temperature and mesh adaptation

Several snapshots of reference case #1 at t = 5 ms when the laser arrives to the position X = 1.2 mm are shown in Figure 3.17. The melt pool is elongated behind the laser spot with a length of about 400 µm. It is not continuous as separate liquid droplets are observed. The formation of droplets will be detailed later. The different iso-contours for temperatures reveal an evolution of the tail shape during cooling. It changes from a convex shape at high temperature to a non-convex shape at low temperature, with an enlargement into the powder bed. At high temperature, for example for T ≈ 1100 • C, the solidified track (position 2) is hotter than its lateral sides (position 1). However, when it is cooled to 700 • C, the temperature in track (position 4) is almost the same as that of its neighboring powder (position 3). This evolution is caused by the insulant effect of powder due to its lower thermal conductivity

(≈ 0.2 W • m -1 • K -1
) as a porous material.

The principle of mesh adaptation was presented in Section 3.4. Here we give more details on its implementation. The first metric is based on the error estimation (Section 3.4.1) with 4 tracked variables: density, liquid fraction, heat source distribution and temperature. By using theses variables, we aim at refining the gas/material interface, the melt pool boundary, the region impacted by laser and the region with high variation of thermal gradient, respectively. The second metric is isotropic in order to refine the region where the powder is melted. The third metric is also isotropic, used to control the global size and avoid the coarsening of elements under a controlled number of total elements. These three metrics are intersected (Section 3.4.2) in a final stage. Figure 3.17(b-d) show the final effect of this mesh adaptation strategy combining 3 metrics. Due to the variation of density gradient (∇{ρ}) across the gas/material interface, the interface is refined and a good representation of the solidified track surface is achieved, as shown in Figure 3.17(b). However, the density field does not provide a refinement to the melt pool boundary as the same density is considered for the solid and liquid. Consequently, we adapt the mesh according to the liquid fraction, giving a refinement of melt pool boundary. Note that at the gas/liquid interface, both the density and liquid fraction make a contribution to mesh refinement. This provides a good tracking of droplet formation (detailed later) as shown in Figure 3.17(c). Different contributions in mesh adaptation are presented in Figure 3.17(d). The refinement in regions A, B and C are due to the density, liquid fraction and heat source fields (first metric), respectively. One can see that the heat source distribution qL respects well the evolution of the gas/material interface, which is indicated by qL = 0 above the interface as α g = 0 for the gas. However, the effect of temperature on mesh adaptation is not obvious as it may be less important than other variables. In the region with high variation of thermal gradient, such as in the powder in front of the laser, the refining effect of temperature field may be superimposed by other contributions. Region D (second metric) is the transition region (exposed to laser irradiation) refined by an isotropic metric of 5 µm mesh size. Its effect can also be seen in (c). The global mesh size is controlled by an isotropic metric of 50 µm as shown in the dense solid material (region E, third metric).

1.2 mm 0.2 X = 0

T l = 2104 • C X Y • 1 • 2 • 3 • 4 (a)

SINGLE TRACK SIMULATION WITH POWDER ON SUBSTRATE

The results in Figure 3.17 are obtained with the blockage of elements (Section 3.4.3) to enforce the representation of the gas/material interface. It is applied around the gas/material interface once the track is solidified (region A) or not impacted by laser (region B). In addition, the powder/substrate interface is also blocked when it is not impacted by laser (region B) in order to avoid coarsening during remeshing. Figure 3.18 shows that the magenta element can be well blocked during mesh adaptation at t = 6, 8 and 10 ms. This is very important and its effect is compared with an additional case without blockage as shown in Figure 3.19 taken at t = 2, 3, and 3.6 ms. During the simulation, the solidified track surface evolves due to the slight evolution without blockage through the successive remeshing and transportation of level set field while it is kept immobile with blockage. Consequently, the blockage is necessary to access the track surface evolution during processing. Note that the number of elements increases with the blockage strategy.

Melt pool dynamics

The droplets that can be observed in Figure 3.17 and 3.19 need further study on their formation. Hence, a series of snapshots is presented in Figure 3.20. When the powder is exposed to the laser irradiation, it is quickly melted due to the high energy concentration of the laser. The resulted liquid is then spheroidized under the effect of surface tension and small droplets are formed (t = 6.12 ms). Note that due to the small dimension (i.e. high curvature) of these droplets, the surface tension is much more significant than the gravity and droplets can be formed. The formed droplets grow as more liquid is formed (t = 6.15 ms). Under the gravity, these droplets fall down into the melt pool, resulting in a liquid wave towards the X Y X = 1.025 mm X = 1.44 mm 6.08 ms 6.12 ms 6.15 ms

SINGLE TRACK SIMULATION WITH POWDER ON SUBSTRATE

tail of the melt pool (t = 6.16 ms). This effect may lead to ripples on the surface of solidified track when the wave arrives to the tail of the melt pool where final solidification occurs [START_REF] Rappaz | Analysis of Ripple Formation in Single Crystal Spot Welds[END_REF].

The formed droplets may furtherly coalesce together to form bigger droplets as indicated at t = 6.24 ms. As shown in this sequence of snapshots, there is a kind of periodicity which can be observed in this phenomenon. It should also be underlined that the fall of droplets due to gravity and the spheroidization due to surface tension generate high velocities (∼0.5 m • s -1 in this example). The solver with semi-implicit formulation of surface tension demonstrates its robustness, delivering stable results to reproduce this dynamics.

Verification of mass and energy conservation

In the simulation with evolving interface tracked by level set and refined by frequent mesh adaptation, the conservation of mass and energy is a challenge. For the conservation of mass, the total mass deviation of the material domain, ∆m, is not representative, neither the ratio between ∆m and the mass of material domain, as it depends largely on the dimension of the substrate. Consequently, we propose to evaluate the ratio between ∆m and the mass of laser impacted powder, m imp , as follows:

∆m m imp = m(t) -m 0 ρ Z 1 2r int ∆h p v L t (3.55)
where m(t) and m 0 are the current and initial mass, respectively. The laser impacted powder volume is approximately estimated by the scanned volume 2r int ∆h p v L t. Note that here the impacted cross section (perpendicular to scanning direction) is calculated by 2r int ∆h p , while the real section differs from this value as it depends on the laser power and scanning speed. However, we aim to give a tendency of mass deviation rather than to evaluate its precise magnitude. The evolution of ∆m/m imp is shown in Figure 3.21, revealing the problem of mass conservation and the normalized deviation tends to be stable after 5 ms. This mass deviation may be related to the remeshing and the calculation of volume contraction rate θ (Eq.3.4), remaining an important point to be optimized in future work.

The conservation of energy for the simple 1D case is verified in Section 3.5 between the theoretical energy E th (calculated with heat input) and the current energy E of the system. Instead of comparing E and E th , here we propose to investigate the energy deviation normalized by the theoretical energy:

∆E E th = E(t) -E th (t) E th (t) (3.56)
where the initial energy of system is set to zero. The result is shown in Figure 3.21(b). At the end of the simulation, E is 1.6% lower than E th , which can be considered admissible.

Influence of material and process parameters

Velocity field

The final track shape is the result of melt pool dynamics combined with solidification. The melt pool dynamics in cases #2-5 is studied by velocity field plotted in Figure 3.22. Case #2 with increased (magnitude) negative Marangoni coefficient ∂γ/∂T shows obvious centrifugal convection flows with high velocity from the hot center to the colder boundary under the laser spot. Two opposite convection cells are formed as shown in the transversal cut view (Figure 3.22(b)). By contrast, the convection cells in case #3 are in different directions than those in case #2, due to the opposite value of ∂γ/∂T , as shown clearly in the zoomed rectangles. Similar to case #2, the generated velocity is higher than that in case #1 (Figure 3.20). Note that case #2 shows a different melt pool shape (transversal cut view) to case #3 near the powder/substrate interface (in red rectangle), which is related to the expansion effect of negative ∂γ/∂T . Another effect of the increasing "negative Marangoni" (case #2) is a longer melt pool than in the reference case #1. Controversely, for a "positive Marangoni", the melt pool is shorter. This is fully consistent with the transport of the high temperature fluid (see also next paragraph). However, in both cases #2 and #3, the Marangoni convection is limited, affecting a small depth due to the rather high liquid viscosity. In case #4 with increased viscosity, a very stable melt pool without droplets is obtained. The resulted velocity is much lower than that in any other cases. The melt pool is always continuous and shows a regular shape. By increasing the scanning velocity in case #5, the linear energy (P L /v L ) is decreased. The shape of melt pool and solidified track shows more irregularity, which is due to the combined effect of the fall of droplets and of solidification. In fact, the fall of droplets in the principal melt pool forms a zone with a larger width, while solidification is so quick that the surface tension has not enough time to smooth the surface. Hence, this larger zone is fixed and necks appear in the solidified track.

SINGLE TRACK SIMULATION WITH POWDER ON SUBSTRATE

expansion of the melt pool width while the length is almost the same. This result may be caused by the combination of two effects. On one hand, the enforced outward (top view) convection flow tends to expand the melt pool. On the other hand, this convection flow provides a better cooling, leading to a smaller melt pool. Consequently, the influence of negative Marangoni coefficient ∂γ/∂T with increased magnitude is not very remarkable in this case. It is more evident with a shorter melt pool [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramic materials[END_REF]. By contrast, with a positive ∂γ/∂T in case #3, both the convection direction and cooling effect result in a smaller melt pool or temperature iso-contours. Hence, their shapes are significantly changed. The melt pool depth is not affected by the modification of ∂γ/∂T as the convection flow is limited to a certain depth as shown previously in Figure 3. 22(b). This may be related to the high viscosity of alumina. One difference is the iso-contour T = 2700 • C, which becomes shallower in case #2 than in case #1. The better cooling in cases #2 and #3 decreases the maximum temperature below to 3100 • C, especially in case #3.

The increased liquid viscosity in case #4 prevents high convection flow. A longer melt pool is observed as the weak convection flow can not efficiently cool the melt pool. The region with temperature higher than 3100 • C is larger than case #1. The melt pool depth is not influenced as the convection flow does not affect the bottom of the melt pool in any cases #1-4 and it is mainly determined by the energy deposition. In case #5 with high scanning speed, due to the decrease of linear energy, a shorter, narrower and shallower melt pool is achieved, which can just melt the substrate.

In order to evaluate the Marangoni effect, the Marangoni number M a is calculated according to Eq.3.15 for cases #1-5 as shown in Table 3.4. According to Maroto et al. [START_REF] Maroto | Introductory analysis of Bénard-Marangoni convection[END_REF], Marangoni convection begins in case #1-3 as |M a| exceeds the critical value 81, while it does not occur (or is too weak) in case #4-5. This distinction qualitatively corresponds to the convection velocity in Figure 3.20 and 3.22. The high M a in case #2 and #3 explains the high velocity obtained in Figure 3

T max • C 3185 3030 2898 3620 2575 ∂γ ∂ T ×10 -5 N • m -1 • K -1 -
= 5.5 W • m -1 • K -1 , ρ = 3970 kg • m -3 and C p = 1900 J • kg -1 • K -1 . Note that the boiling temperature of Al 2 O 3 is 2967 • C.

Track surface morphology

The track surface morphology is the direct output of combined melt pool dynamics and solidification. High melt pool dynamics results in irregular track surface, even with interuption of track or balling effect. The track surface morphology can be visualized by height isocontours as shown in Figure 3.24. Note that the dimension in the horizontal Y direction is scaled by a factor 2 in order to better present the iso-contours, while there is no scale change in X and Z directions. At the beginning of the track, a transition region is observed. However, the domain of interest is the quasi-steady state which follows. Firstly, it should be recognized that continuous iso-contours mean smooth track surface, while it has more surface roughness when iso-contours are separated as indicated for case #1. At the same time, high gradient of color level also indicates high surface roughness. Hence, case #2 suggests less surface roughness than case #1, while opposite effect is observed in case #3. This can be explained by the opposite Marangoni convection in cases #2 and #3. As the regions with high altitude (hills) is usually caused by the fall of hot droplets, negative Marangoni coefficient ∂γ/∂T moves the liquid from hills to valleys (low region), thus smoothing the surface. Consequently, a positive ∂γ/∂T intensifies this height difference. The melt pool dynamics decreases with higher viscosity in case #4, leading to a much smoother track surface. In case #5, the affecting duration of surface tension is decreased by two as the solidification velocity (∼ v L ) is doubled. As a result, the melt pool can not be smoothed under surface tension before solidification and high surface roughness is generated.

CONCLUSION

Conclusion

Under the assumption of a continuous powder bed, the numerical modeling of melt pool dynamics is established by coupling mass and momentum conservation, with a compressible Newtonian behavior for the powder and an incompressible Newtonian behavior for liquid and air. The surface evolution of the melt pool can be tracked by the level set surface ψ = 0. The modeling and the implementation of driving forces are essential to the melt pool dynamics. Among them, constant gravity, surface tension and Marangoni forces are taken into account.

Although other forces like recoil pressure are also important, they are not yet integrated in the present model. The FE implementation is presented, particularly in the semi-implicit formulation of surface tension. Mesh adaptation is discussed as it is important in the model with significant surface evolution.

The model is firstly validated by a simple case. Then the influence of surface tension, viscosity and scanning speed on the melt pool dynamics is investigated by single track simulation with powder deposition on a substrate. The Marangoni effect can eventually change the convection flow inside the melt pool and its shape, with different influences depending on the sign of the Marangoni coefficient ∂γ/∂T . This effect can redistribute heat and decrease the melt pool temperature. In addition, the sign of ∂γ/∂T has a direct influence on the track surface morphology. Another important material property is the liquid viscosity, which has an opposite effect than to surface tension. The higher is the viscosity, the less surface roughness is obtained. Calculations show that the developed simulation is able to detect the occurrence of balling effect, where the track surface quality is degenerated when the speed attains certain limit.

Résumé

La mécanique des solides est modélisée par un modèle élasto-viscoplastique. Le lit de poudre est assimilé à un fluide newtonien car on ne s'intéresse pas à la mécanique dans le lit de poudre. L'avantage de ce modèle est sa généralité. Il peut être facilement appliqué à d'autres matériaux. L'implémentation éléments finis et la résolution numérique sont ensuite détaillées. Le changement de variable sur la variable pression est présenté dans le contexte d'une résolution mixte vitesse/pression. Son optimisation significative sur la réduction du temps de calcul est soulignée.

Les simulations sont présentées pour étudier l'effet du laser auxiliaire sur la distribution des contraintes. Des contraintes de traction élevées dans la direction de balayage (x) peuvent être générées dans la zone refondue derrière le bain liquide, conduisant à des fissures transversales potentielles. Les contraintes de traction dans les directions transversale (y) et de construction (z) sont responsables des fissures potentielles dans le plan longitudinal et le délaminage, respectivement. Elles sont beaucoup plus faibles que dans la direction de balayage. De plus, une rugosité de surface élevée peut créer une concentration des contraintes, ce qui augmente le risque de fissuration. L'utilisation d'un laser auxiliaire peut réduire la contrainte maximale. Un laser auxiliaire à l'arrière du bain liquide est plus efficace pour diminuer la contrainte maximale, en réduisant le gradient thermique dans la zone de refroidissement. Le type de source de laser auxiliaire est également important. Le laser auxiliaire avec une pénétration plus profonde (Yb:YAG) est meilleur que le laser limité à une action de surface (CO 2 ) pour réduire la contrainte maximale. Le code de calcul développé pourra être utilisés pour choisir un laser auxiliaire et une stratégie de chauffage afin de produire des pièces sans fissures.

As thermally induced stresses and strains in SLM process are responsible for cracks and distortion, solid mechanics will be modeled in this chapter during track formation. After a literature review on the modeling of solid mechanics, the model used in this work will be presented, particularly a general elasto-viscoplastic model in the context of multiphase domain with level set method. The FE implementation and numerical resolution are then detailed, together with a technique of computing acceleration. The model is used to simulate the stresses in single track development. The investigation focuses on the influence of an auxiliary laser on the maximum stress during processing, including the type of auxiliary laser and its heating strategy.

State of the art

Modeling of solid mechanics in additive manufacturing provides access to the evolution of stress and strain distributions during processing. This functionality is available in commercial softwares like Abaqus [START_REF] Vastola | Controlling of residual stress in additive manufacturing of Ti6Al4V by finite element modeling[END_REF], Ansys [START_REF] Zhang | Finite element simulation and experimental validation of distortion and cracking failure phenomena in direct metal laser sintering fabricated component[END_REF], Simufact Additive, Comsol and SYSWELD [START_REF] Wang | Residual stresses in LENS-deposited AISI 410 stainless steel plates[END_REF] or the AM package of ESI group. Finite element method is usually used in commercial software and often accompanied with element activation or deactivation like in Comsol. Most commercial software pay more efforts to the computing efficiency and the modeling at part scale. Many researches [START_REF] Zhang | Finite element simulation and experimental validation of distortion and cracking failure phenomena in direct metal laser sintering fabricated component[END_REF][START_REF] Kamara | A Numerical Investigation into Residual Stress Characteristics in Laser Deposited Multiple Layer Waspaloy Parts[END_REF][START_REF] Megahed | Metal additivemanufacturing process and residual stress modeling[END_REF] have been conducted at part scale to predict distortion and residual stress. With the thermomechanical model, Zhang et al. [START_REF] Zhang | Macro-and mesoscale finite element modelling of the laser beam melting process[END_REF] predict the stress distribution during construction. Their results in Figure 4.1 show that σ xx (in the construction plane) is higher than σ zz (in the building direction) and high σ xx can be generated in the substrate. However, simulations at part scale are often carried out with predefined track shape as with a flat surface, leading to the loss of precision of stress prediction.This type of investigation can help to optimize the part geometry and process parameters. On the other hand, the modeling of solid mechanics at track scale can provide more precise stress distribution and evolution around melt pool during heating and cooling. Hence, potential cracks can be predicted and heating strategy can be adapted to avoid the occurrence of this default. In the following, the literature review will focus on modeling at track scale. Gusarov et al. [START_REF] Gusarov | Residual Stresses at Laser Surface Remelting and Additive Manufacturing[END_REF] studied the stress distribution with a thermoelastic model proposed for laser surface remelting or SLM processing. They generalized the Hooke's law by introducing an additional isotropic tension term inside the remelted zone. This term represents the complete relaxation of the compressive stresses at temperatures above the melting point. Some important features of this model are as follows:

• viscoplasticity is absent; • transverse section shape of remelted zone is assumed to be known a priori and uniform in the scanning direction (x) as indicated by the schematic view of simulated system in They investigated the influence of remelted zone on stress distribution with SiO 2 . The influence of initial temperature was also studied with Al 2 O 3 . The residual stresses are found to be 4.1. STATE OF THE ART dependent on the shape of remelted zone but independent on its size. The maximum tensile stress in scanning direction is about twice than that in transverse direction, while it is much smaller in the building direction. These findings support the principal crack formation observed in experiments as shown in Figure 1.8. In addition, the maximum tensile stress is not influenced by the shape of the remelted zone. Another useful finding is the significant effect of preheating temperature on the reduction of maximum stress. Inspired by Hagedorn [START_REF] Hagedorn | Net shaped high performance oxide ceramic parts by selective laser melting[END_REF] et al., the preheating of Al 2 O 3 is simply approached by the choice of initial preheating temperature T 0 equal to 1600 • C in the simulation. The result is compared with T 0 = 0 • C as shown in Figure 4.2. The maximum tensile stress occurs at the boundary of remelted zone. The preheating can largely decrease the maximum value by a factor 5 in the scanning direction x (from 7.4 to 1.7 GPa) and in the transverse y (from 5 to 1 GPa) and build direction z (from 0.3 to 0.06 GPa). Consequently, this preheating appears to be very useful for eliminating cracks in SLM process. Although this model predicts qualitative stress distributions, it may be too simple to provide precise results due to its assumption listed above. Particularly, the shape of the remelted zone is certainly not uniform in real process and our simulations will show that it has important effect on stress distribution. Vastola et al. [START_REF] Vastola | Controlling of residual stress in additive manufacturing of Ti6Al4V by finite element modeling[END_REF] studied the influence of laser size, power, scanning speed and initial temperature on the stress distribution in single track processing of Ti -6Al -4V with EBM. They used FE method in Abaqus commercial software taking into account elasto-plastic and thermal strains. A multiphase model was implemented, including powder, liquid and solid, with temperature dependent material properties. Among them, powder and liquid were modeled as perfectly plastic materials with negligible yield strength for the sake of better numerical convergence. By solving the heat transfer equation, temperature field was obtained and used to deduce the thermal expansion or shrinkage. The stress distribution was predicted under the assumption of flat melt pool surface. The results show that, regardless of the influence of laser size, power and scanning speed on the depth of remelted zone, they have no significant effect on the maximum tensile stress located in the remelted zone as shown in Figure 4. 3(a-c). However, each increase of initial temperature by 50 • C leads to a stress reduction of ∼ 20% according to Figure 4.3(d).

With the modeling of heat transfer and melt pool dynamics presented in previous chapters, we can model the solid mechanics with temperature distribution and shape of melted zone. Hence, more details can be accessed like the influence of track height on stress distribution. On the other hand, the literature results lead us to focus the investigation on the influence of the preheating strategy on stress distribution, rather than the influence of laser parameters.

CHAPTER 4. SOLID MECHANICS

Depth from free surface ( µm)

σ xx (MPa) (a)
Depth from free surface ( µm)

σ xx (MPa) (b)
Depth from free surface ( µm)

σ xx (MPa) (c)
Depth from free surface ( µm)

σ xx (MPa) (d)

Modeling of solid mechanics

In SLM process, the stresses are generated by thermal expansion or shrinkage due to the density variation during heating or cooling, respectively. The expansion or shrinkage has thermal effect in a point of view of energy conservation. However, this effect is negligible compared to the laser energy input [START_REF] Gusarov | Thermoelastic Residual Stresses and Deformations at Laser Treatment[END_REF]. Therefore, in the numerical model presented hereafter, the solid mechanics is not strongly coupled with thermal modeling. This numerical model is principally based on previous works at Cemef by Jaouen [START_REF] Jaouen | Modélisation Trdimensionnelle Par Éléments Finis Pour l'analyse Thermo-Mécanique Du Refroidissement Des Pièces Coulées[END_REF], Hamide [2] and Desmaison [START_REF] Desmaison | Modélisation Numérique d'un Procédé de Soudage Hybride Arc / Laser En Approche Level Set : Application Au Soudage Multi-Passes de Tôles d'acier de Forte Épaisseur[END_REF]. In fact, modeling of solid mechanics in SLM is similar to the one developed for welding process presented by Hamide and Desmaison despite the presence of powder in SLM.

Governing equations

Similar to the modeling of fluid dynamics in the previous chapter, the modeling of solid mechanics consists of coupling the conservation equations of momentum and mass. In the context of solid mechanics, we neglect in a first approach the inertia effect, leading to the equilibrium equation:

-∇ • σ = f v (4.1)
where f v = ρg. The boundary conditions are:

     u = u imp on Γ u σ • n Ω = f imp on ∂Ω f (4.2)
It should be mentioned that the resolution of equilibrium equation is mainly governed by the constitutive law describing the relation between stress and strain.

Constitutive law

In the context of additive manufacturing or welding, material endures a large temperature variation from ambient temperature up to melting point. The modeling of solid mechanics requires a suitable constitutive law which can be used in the full temperature range. As viscous phenomena become important at high temperature, a general elasto-viscoplastic model is necessary [2]. In the present model, 3 regimes are distinguished:

• elasto-viscoplastic (EVP) behavior at temperature T lower than a critical temperature T c (usually chosen close to solidus T s ) when viscous effects are still negligible; • viscoplastic (VP) behavior in the mushy zone when T c ≤ T ≤ T l (liquidus); • newtonian behavior when T > T l .

Similar to the model precedently used in welding by Hamide [2] and Desmaison [START_REF] Desmaison | Modélisation Numérique d'un Procédé de Soudage Hybride Arc / Laser En Approche Level Set : Application Au Soudage Multi-Passes de Tôles d'acier de Forte Épaisseur[END_REF], the first regime can be modeled by an EVP constitutive law while other two behaviors can be both described by a VP constitutive law as the Newtonian behavior is just a special case of VP behavior.

In the case of small deformation, the strain rate ˙ can be decomposed into several contributions:

˙ =      ˙ th + ˙ el + ˙ vp if T < T c ˙ th + ˙ vp if T > T c (4.3)
where ˙ th , ˙ el , ˙ vp are the thermal, elastic and viscoplastic parts of ˙ , respectively. The thermal expansion rate ˙ th can be directly related to the density variation by:

˙ th = - 1 3ρ dρ dt I (4.4)
Taking ė as the deviatoric part of ˙ , ˙ can be decomposed into deviatoric and spherical parts (example for EVP):

     ė = ėel + ėvp tr ˙ = tr ˙ el + tr ˙ th (4.5)
as ˙ th is purely spherical and tr ˙ vp = 0. Similarly, σ is also decomposed as follows:

σ = s -pI (4.6)
where p = -1 3 tr σ is the pressure. The stress-strain relation for each constitutive law will be detailed in the following.

Elasto-viscoplasticity

For the elasto-viscoplastic constitutive law, the elastic strain el can be linearly related to stress σ by the Hooke's law:

σ = E : el (4.7)
where E is the 4th order symmetric stiffness tensor with 21 independent components [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF].

The elastic strain rate ˙ el can be then obtained:

˙ el = E -1 σ + ∂E -1 ∂t σ (4.8)
In the case of isotropic material, E depends only on Young's modulus E and Poisson's ratio ν, so Eq.4.8 can be written as:

˙ el = 1 + ν E σ - ν E tr σ I + d dt 1 + ν E σ - d dt ν E tr σ I (4.9)
By splitting the spherical and deviatoric part, the expression of the elastic part of the strain rate tensor is:

     ėel = ṡ 2µ -1 2µ 2 dµ dt s tr ˙ el = -ṗ χ b + 1 χ 2 b dχ b dt p (4.10) 102

MODELING OF SOLID MECHANICS

where µ = E 2(1+ν) is the shear modulus and χ b = E 3(1-2ν) is the bulk modulus. A multiplicative law is used for the expression of the yield stress. The equivalent stress is then written as:

σ = σ Y + K √ 3 m+1 ˙ m ¯ n (4.11)
where σ Y is the nominal yield stress, K the consistence coefficient, m the sensibility coefficient and n the hardening coefficient. They are all assumed to be only temperature dependent. σ = 3 2 s : s is the equivalent stress of von Mises and ˙ = 2 3 ˙ vp : ˙ vp is the equivalent viscoplastic strain rate. The relation between strain rate and the deviatoric part of stress can be then obtained:

˙ vp = 3 2 ˙ σ s = √ 3 2σ σ -σ Y √ 3K¯ n 1/m + s = λs (4.
12)

The function f + takes the value of f when f is positive and takes 0 otherwise. This means that viscoplasticity is evoked only when the equivalent stress σ exceeds the yield stress. The accumulated plastic strain ¯ can then be calculated by:

¯ = t 0 ˙ dt (4.13)
Combining Eq.4.4, 4.5, 4.12 and inserting them into Eq.4.10, the EVP constitutive law can be written as:

       ṡ = 2µ ė -λs + 1 µ dµ dt s ṗ = -χ b tr ˙ + 1 ρ dρ dt + 1 χ b dχ b dt p (4.14)

Viscoplasticity

The Norton-Hoff model is used for the purely viscoplastic constitutive law:

˙ vp = 1 2K √ 3 ˙ 1-m s (4.15)
leading to the yield stress expressed as:

σ = K √ 3 m+1 ˙ m (4.16)
Note that Eqs.4.15 and 4.12 model have the same expression when σ Y → 0, m → 1 and n → 0. In this particular situation, Eq.4.15 is simplified to:

˙ vp = 1 2K s (4.17)
Thus the Newtonian behavior is retrieved with a consistence coefficient K equivalent to the dynamic viscosity η.

Numerical resolution

Eq.4.1 coupled with material constitutive law is solved by finite element method, particularly with a P 1 + /P 1 formulation. This is implemented in Cimlib and the corresponding solver named TransWeldEVPMini. As the stress-strain relation is nonlinear, the search of equilibrium state of whole system uses the Newton-Raphson iteration. In order to reduce the computation time of each iteration, the solver gets benefits from the change of variable.

When calculating the contribution of each element in stiffness matrix during each iteration, the stress-strain relation is determined by one of the constitutive laws detailed in the previous section. In the case of EVP constitutive law, this also needs local resolution of a nonlinear scalar equation with Newton-Raphson iteration to find the unknown λ, or equivalently ˙ . All these numerical treatments will be detailed in the following, together with an algorithm flow chart. At the end, a global flow chart will be presented to explain the coupling of heat transfer, fluid mechanics and solid mechanics.

FE implementation

Inserting Eq.4.6 and the second equation of Eq.4.14 into Eq.4.1 leads the equilibrium and mass conservation equations:

     ∇ • s -∇p + f v = 0 ∇ • u + ṗ χ b -1 χ 2 b dχ b dt p + 1 ρ dρ dt = 0 (4.18)
As mentioned in the previous chapter for the FE implementation of momentum and mass conservation equation, the instability with simple P 1 formulation for both velocity and pressure can be overcome by the P 1 + /P 1 formulation with additional degree of freedom introduced by bubble function in velocity field. After the discretization of domain Ω by

Ω h = E∈T h Ω E
(E is any element of mesh T h ), the following functional spaces are then defined:

V h = u h |u h ∈ P 1 Ω E d , u h = u imp on Γ h u (4.19) V h0 = u h0 |u h0 ∈ P 1 Ω E d , u h0 = 0 on Γ h u (4.20) B h = b h |b h ∈ P 1 Ω E d , b h = 0 on ∂Ω E (4.

21)

P h = p h |p h ∈ P 1 Ω E (4.22)
The approximated velocity and pressure fields are:

w h (x) = u h (x) + b h (x) = N i (x)U i + N b (x)B (4.23) p h (x) = N i (x)P i (4.24)
where N i (x) (i ∈ [1, D] the local numbering of node in an element) is the linear interpolation function and N b (x) the bubble function. Here Einstein's notation is applied with dummy indices to be summed. Note that U i and N i are associated to each node i while B and N b are associated to element E. An element E can be decomposed to D sub-elements E i connected to barycenter G as illustrated in

i j k G E i E j E k (a) i j k G N b (x) 1 ξ η (b) G j k 1 i ξ η N i (x) N b (x) (c)
Ω E ∇N b (x) dΩ = Ω E i D ∂N i (x) ∂ x j dΩ = D V E D i ∂N i (x) ∂ x j = 0 (4.25)
as ∂N i (x)/∂x j is constant in each sub-element and

i N i (x) = 1.
Here V E is the volume of element E. This means that for any constant tensor T :

Ω E T : ∇b h dΩ = T : B ⊗ Ω E ∇N b (x) dΩ = 0 (4.26)
The FE implementation corresponds to find (

w h = u h + b h , p h ) ∈ (V h ⊕ B h , P h ), so that for ∀ w * h = u * h + b * h , q * h ∈ (V h0 ⊕ B h , P h ):            s(w h ) : ˙ w * h -p h , ∇ • w * h -f imp , w * h Γ h f -f v , w * h = 0 -∇ • w h , q * h + 1 χ 2 b dχ b dt -1 χ b ∆t p h , q * h - 1 ρ dρ dt - p - h χ b ∆t , q * h = 0 (4.
                   s(u h + b h ) : ˙ (u * h ) -p h , ∇ • u * h -f imp , u * h Γ h f -f v , u * h = 0 s(u h + b h ) : ˙ (b * h ) -(p h , ∇ • b * h ) -f imp , b * h Γ h f -(f v .b * h ) = 0 -∇ • (u h + b h ), q * h + 1 χ 2 b dχ b dt -1 χ b ∆t p h , q * h 1 ρ dρ dt - p - h χ b ∆t , q * h = 0 (4.28)
The red term in Eq.4.28 vanishes as the bubble function vanishes at the boundary of each element. Jaouen [START_REF] Jaouen | Modélisation Trdimensionnelle Par Éléments Finis Pour l'analyse Thermo-Mécanique Du Refroidissement Des Pièces Coulées[END_REF] showed that s(u h + b h ) can be decomposed into s(u h ) + s(b h ). Hence, with some arrangement, Eq.4.28 can be formatted to:

                 s(u h ) : ˙ (u * h ) + s(b h ) : ˙ (u * h ) -p h , ∇ • u * h -f imp , u * h Γ h f -f v , u * h = 0 s(u h ) : ˙ (b * h ) + s(b h ) : ˙ (b * h ) -(p h , ∇ • b * h ) -(f v , b * h ) = 0 -∇ • u h , q * h -∇ • b h , q * h + 1 χ 2 b dχ b dt -1 χ b ∆t p h , q * h -1 ρ dρ dt - p - h χ b ∆t , q * h = 0 (4.29)
Note that the red terms in Eq.4.29 vanish according to the property of bubble function in Eq.4.26. This is the system to be solved and it can be formulated in a residual form:

           R U = R uu + 0 + R up -F u = 0 R B = 0 + R bb + R bp -F b = 0 R P = R pu + R pb + R pp -F p = 0 (4.30)
with corresponding terms in the same color as in Eq.4.29.

Non-linear resolution

The system in Eq.4.30 is non-linear. Newton-Raphson method is proposed for its resolution, which means that we should firstly obtain the stiffness matrix K with components calculated by:

K xy = ∂R x ∂ y x, y ∈ {U , B, P } (4.31)
This results in the final system to be solved:

      K U U 0 K U P 0 K BB K BP K P U K P B K P P             δU δB δP       = -       R U R B R P       (4.32) 106

NUMERICAL RESOLUTION

As we are not interested in the additional velocity field B and due to the property of bubble function, it can be eliminated by condensation in each element (see more details in [START_REF] Jaouen | Modélisation Trdimensionnelle Par Éléments Finis Pour l'analyse Thermo-Mécanique Du Refroidissement Des Pièces Coulées[END_REF]):

δB = -K BB -1 R B + K BP δP (4.33)
Inserting Eq.4.33 into Eq.4.32, the bubble function can then be eliminated:

   K U U K U P K P U K P P -C       δU δP    = -    R U R pu + R pp + F p -CP t -K P B K BB -1 F b    (4.34)
with C = K P B K BB -1 K BP . The solution of velocity and pressure variations is then used to update the solution from the previous iteration (υ -1) to the current iteration (υ):

   U P    (υ) =    U P    (υ-1) +    δU δP    (4.35)
The computing of the stiffness matrix requires the derivative of deviatoric stress with respect to the velocity:

∂s ∂ U = ∂s ∂ ˙ ∂ ˙ ∂ U (4.36)
The second part at RHS can be calculated based on the linear interpolation function [START_REF] Rappaz | Numerical Modeling in Materials Science and Engineering[END_REF]:

∂ ˙ ∂ U ijkl = 1 2 ∂N k ∂ x j δ il + ∂N k ∂ x i δ jl (4.37) with i, j, l ∈ [1, d] and k ∈ [1, D].
The first part ∂s/∂ ˙ is called tangent modulus and its calculation should be addressed to the constitutive law. This will be discussed in the next section.

Local resolution

The stress-strain state in each element is the unknown to be find and also necessary for the calculation of tangent modulus. Based on the Eq.4.14 and using a fully implicit schema, we can get:

s t -s t-∆t ∆t = 2µ t ėt -λt s t + µ t -µ t-∆t µ t-∆t ∆t s t-∆t (4.38)
leading to:

s t = µ t µ t-∆t s t-∆t + 2µ t ėt ∆t 1 + 2µ t λt ∆t (4.39)
Here all notation at current time step t correspond also to the current iteration (ν) by default. Note that if λt = 0 (pure elastic), Eq.4.39 gives the elastic prediction of s t . Another relation can be established from :

s t : s t = 2 3 σ2 (4.40)
Inserting Eq.4.39 into Eq.4.40 and using Eq.4.12, we can finally get:

3µ t ˙ t ∆t + σ(¯ t-∆t , ˙ t ) -B 0 = 0 (4.41)
where

B 2 0 = 3 2
µ t µ t-∆t s t-∆t + 2µ t ėt ∆t :

µ t µ t-∆t s t-∆t + 2µ t ėt ∆t (4.42)
Eq.4.41 is a non linear equation and it is solved by Newton-Raphson iteration to find ˙ t . Then ¯ t can be obtained by Eq.4.13. With Eq.4.11, the equivalent stress σ can be calculated. Later we can get λt by the flow law in Eq.4.12 and finally get s t by Eq.4.39.

The stress-strain state is fully deduced at this stage, making it possible to calculate the tangent modulus for each element:

C = ∂s t ∂ ˙ t =   2µ σ B 0 I -I ⊗ I 1 d -2µs ⊗ s 3 2σ 2   σ B 0 -1 + 1 1 + ∂ σ ∂¯ 1 3µ      ∆t (4.43)
or we can use the index notation for tensors:

C = C mnpq I = δ mp δ nq I ⊗ I = δ mn δ pq s ⊗ s = s mn s pq (4.44)
where δ mp is the Kronecker delta. The detailed demonstration can be found in the work of Gay [START_REF] Gay | Contribution à La Simulation Numérique Tridimensionnelle Du Forgeage à Froid[END_REF]. The tangent modulus related to the bubble function is:

C b = 2µ σ B 0 I -I ⊗ I 1 d ∆t (4.45)
Note that both the stress-strain state and tangent modulus are deduced in each element before the construction of local matrix. They appear in K U U and K BB in Eq.4.34.

Change of variable

Eq.4.34 can be simplified by noting:

   K (U U ) K (U P ) K (P U ) K (P P )    (υ)    δU δP    = -    R (U ) U (υ-1) , P (υ-1) R (P ) U (υ-1) , P (υ-1)    (4.46)
In this mixed resolution for velocity/pressure, the diagonal elements of K (U U ) and K (P P ) have very different orders of magnitude. This leads to high condition number, said illconditioned, which has severe consequence in the convergence of linear resolution.

The change of variable was employed by Gay [START_REF] Gay | Contribution à La Simulation Numérique Tridimensionnelle Du Forgeage à Froid[END_REF] to overcome this problem in the elastoplastic modeling. It was later adapted by Jaouen [START_REF] Jaouen | Modélisation Trdimensionnelle Par Éléments Finis Pour l'analyse Thermo-Mécanique Du Refroidissement Des Pièces Coulées[END_REF] in the VP/EVP hybrid modeling in

NUMERICAL RESOLUTION

Thercast. This method consists of changing the pressure p to q: p = q (4.47) leading Eq.4.46 to:

   K (U U ) K (U Q) K (QU ) K (QQ)    (υ)    δU δQ    = -    R (U ) U (υ-1) , Q (υ-1) R (Q) U (υ-1) , Q (υ-1)    K (U Q) = K (QU ) = K (U P ) K (QQ) = 2 K (P P ) δQ = δP / R (Q) = R (P ) (4.48)
The choice for the value of is delicate in order to increase the resolution speed. An optimized value of was proposed by Bellet [START_REF] Bellet | Optimisation du changement de variable sur la pression pour les formulations éléments finis vitesse-pression de THERCAST et R2SOL[END_REF] based on the diagonal components of the elementary stiffness matrix associated to node i:

i = max j∈[1,d] K (U U ) ij,ij K (P P ) ii (4 .49)
Taking i = 1 for all nodes i at initial state, the system is later assembled with new variable Q. Hence, the calculation of should be based on Q. Inserting K (QQ) = 2 K (P P ) into Eq.4.49, at t can be taken as:

t i = t-∆t i max j∈[1,d] K (U U ) ij,ij K (QQ) ii (4.50)
The use of can bring the diagonal components of K (U U ) and K (QQ) to the same order of magnitude. Consequently, the conditioning of the stiffness matrix is improved and the computation time of linear resolution can be significantly reduced. A comparison will be given later regarding the number of iterations and computation time of linear resolution.

Resolution algorithm for solid mechanics

Figure 4.5 shows the algorithm flow chart for the resolution of solid mechanics. At each time step t, the contribution of each element in the global stiffness matrix K and residual R is evaluated by looping over every element. For each element, the solver TransWeldEVPMini is called, corresponding to the flow chart at right in the blue rectangle. The contribution of each element is then assembled into K and R. The coefficient is calculated if it is the first iteration and saved to be used in the next time step. This means that the same (calculated from the previous time step) is used in the current time step. The solution δU and δQ are obtained by solving Eq.4.48 and converted to δU and δP . The velocity and pressure fields are then corrected by this variation. Finally the convergence condition is checked to decide the termination of this time step. Eq.4.50

(U P ) (υ-) , s t-, ¯ t- ˙ , ė T > T c T ⇒ E, ν, µ, χ b , K, m, n λ = 0 ⇒ s elP Eq.
K (U U ) , K (U Q) , K (QQ) R (U ) , R (Q) End N Y N Y N Y Y N TransWeldEVPMini Figure 4
.5: Algorithm flow chart of solid mechanic solver. For the sake of simplicity, the previous time step and iteration are simplified to t-and (ν-), respectively. If not specified, all notations are at current time step t and current iteration (ν).

NUMERICAL RESOLUTION

The stress-strain state in each element is updated during the calculation of elementary contribution. With s t-, ¯ t-at previous time step and (U P ) (υ-) as input, the strain rate ˙ and its deviatoric part ė can be evaluated. According to the local temperature at the barycenter of element, EVP or VP constitutive law can be applied to the element. For example, in the case of EVP behavior, the elastic prediction will be firstly tried by assuming λ = 0 in Eq.4.39. The coefficient B 0 (Eq.4.42) is then compared with σ (Eq.4.11) to check if the element is in pure elastic or elasto-viscoplastic state. If it is elasto-viscoplastic, the equivalent viscoplastic strain rate ˙ is calculated according to Eq.4.41 with Newton-Raphson method. σ, λ and s can then be obtained by Eq.4.11, 4.12 and 4.39. The tangent modulus is finally calculated by Eq.4.43 and 4.45 in order to evaluate the components of stiffness matrix. Firstly, it should be pointed out that at each time step, all the heat transfer, fluid dynamics, transport and reinitialization of LS, solid mechanics and update of coordinates are solved successively. Remeshing is executed with a fixed frequency but not at each time step. The frequency is chosen in the way that the gas/material interface does not exceed the refined transition zone in order to assure good interface representation. Thus, it depends on the time step, the velocity of interface and the size of refined transition zone. Note that it may be adapted according to different melt pool dynamics, for example, when processing metals with much lower viscosity.

Global resolution algorithm

The heat source qL should be firstly obtained as an input for heat transfer. This part is detailed in Section 2.5.1. Together with the convection velocity u - f and temperature T - of previous time step as input, the current temperature distribution can be obtained via the temperature solver (Section 2.5.2). Phase fractions and material properties are then calculated according to temperature. They are used to evaluate the volume variation rate θ (Eq.3.5), Marangoni force f m (Eq.3.13), total curvature κ t (Eq.3.10) and normal direction n (Eq.3.11) related to surface tension force (semi-implicitly calculated inside the solver). As the resolution of fluid velocity is carried out in the whole system, the solid viscosity (1000 Pa • s) is chosen to be just high enough to fix the solid while not too high in order to deliver good resolution convergence. The resulting fluid velocity u f is used to update the level set field by transport and reinitiallization (Section 2.2.1). Particular attention should be paid to the non-conservation character of the reinitiallization method discussed in Section 3.4.3. The resolution of solid mechanics is just after the transport and reinitiallization. Contrary to the fluid mechanics, liquid viscosity will be artificially increased (see the third simplification in Section 4.5) to accommodate the numerical resolution. The results are stress-strain state and velocity, noted as u s in order to distinguish it from u f resulting from fluid mechanics. u s serves to update mesh node coordinates x while keeping its connection T h . At the end, if necessary, the mesh adaptation is executed, including the blockage of elements (Section 3.4.1). This will change the mesh node coordinates and its connection at the same time. Note that in this sequential algorithm, the resolution of heat transfer, fluid dynamics and solid mechanics are weakly coupled, meaning that there is no global iteration procedure over all solvers in each time step.

Initialization

Heat transfer

Fluid dynamics

LS transport + reinitialization

Solid mechanics

Update coordinates

Remesh

t = t end End qL , T -, u - f u - f , κ t , n, f m , θ u f , ψ - T, T -, u - s , p - s , s -, ¯ - u s , x - x -, T - h , M T u f , p f ψ u s , p s , σ, x x, T h t ← t + ∆t Y N

Application to static welding

In order to investigate the effect of change of variable (Section 4.3.4), a benchmark case previously validated by Desmaison [3] is reused to compare the linear resolution time with and without change of variable. This case simulates the static TIG fusion of a metal disk with 30 mm of radius and 10 mm of thickness (Figure 4.7). Material used is the 18M5 steel with thermal and mechanical properties presented in the work of Costes [START_REF] Costes | Modélisation Thermomécanique Tridimensionnelle Par Éléments Finis de La Coulée Continue d'aciers[END_REF]. A surface heat source is imposed at the center of disk during 10 s with a Gaussian distribution:

qL = P L 2πr 2 L exp - r 2 2r 2 L ( 4.51)
with P L = 1200 W and r L = 3 mm. A level set between gas and metal is used with a transition zone of 0.5 mm. The mesh is illustrated in The difference between these two cases is the computing performance as reported in Table 4.1 with the number of iterations (resolution of linear system) for the first Newton-Raphson iteration (υ = 1) at different time steps and corresponding computation time. In fact, the case without change of variable is equivalent to set i = 1. For the first time step, as i is set to be 1 for the case with change of variable, the number of iterations and computation time are almost the same in both cases (first row in Table 4.1). However, the change of variable begins to work from the second time step. It can largely reduce the number of iterations and thus the computation time. The decrease ratio varies with times steps and the final total time of linear resolution is decreased by a factor of about 30. This result shows clearly the benefit of using the change of variable.

Stress distribution in SLM

In this section, the presented model for solid mechanics will be used to simulate the stress distribution during SLM process with alumina. The mechanical behavior of powder is not of interest, neither the gas. Consequently, the constitutive law for different materials in their possibly different states is treated as follows:

• powder and gas are both considered identically as Newtonian fluids with the same value of their viscosity equal to that of liquid alumina; • EVP or VP behavior is applied only when the fraction of dense matter g Z 2 (at element barycenter) attains 1, otherwise Newtonian behavior is applied; • critical temperature T c is set to be 10 • C below the solidus

T s = 2004 • C.
In order to reduce the numerical difficulty, some simplifications proposed by Desmaison [START_REF] Desmaison | Modélisation Numérique d'un Procédé de Soudage Hybride Arc / Laser En Approche Level Set : Application Au Soudage Multi-Passes de Tôles d'acier de Forte Épaisseur[END_REF] are also taken hereafter:

• the density is assumed to be constant when T > T l ; • the sensibility coefficient is taken to be 1 once material attains the critical temperature, meaning that Newtonian behavior is directly present after EVP behavior; • the consistence coefficient (or viscosity) of liquid alumina, gas and powder is increased to 1000 Pa • s; • coordinate update by the velocity resulting from the solution of the solid mechanics is deactivated due the negative element volume caused by element distortion (can be confronted by increasing the liquid viscosity, see Section 5.4). Although Al 2 O 3 is brittle and the viscoplastic behavior seems not appropriate, this general model is employed in the following simulations. It requires mechanical properties, particularly the density ρ(T ), yield stress σ Y (T ), Young's modulus E(T ), Poisson's ratio ν(T ), consistence coefficient K(T ), sensibility coefficient m(T ) and hardening coefficient n(T ). The density of alumina evolves from 3970 to 3780 kg • m -3 with temperature varying from 20 to 2104 • C [START_REF] Aksay | [END_REF]. Sánchez-González et al. [START_REF] Sánchez-González | Application of Hertzian Tests to Measure Stress-Strain Characteristics of Ceramics at Elevated Temperatures[END_REF] characterized the stress-strain relation of alumina from 25 to 1200 • C using Hertzian test. They used a Lüdwig constitutive law to model the stress-strain relation. This is an elasto-plastic model taking into account hardening. By assuming m = 0, we can interpolate their model by the general EVP model and get the stress-strain curve as shown in Figure 4.9(a). In order to introduce a viscous behavior, the evolution of m with respect to the critical temperature T c is taken as the same as in the static TIG welding (Figure 4.9(b)). It should be mentioned that the calibrated model is based on compression test while the tensile strength of alumina is usually between 260 and 300 MPa at room temperature, much lower than the compression strength. Consequently, non-physical tensile stress higher than tensile strength may be predicted in simulation. In fact, it can be considered as a signal of cracks.

In the following, we aim at investigating the effect of additional heating on the stress distribution by several cases listed in Table 4.2. The process parameters of the principal laser are: P L = 84 W, r int = 50 µm and v L = 0.3 m • s -1 . The additional heating is realized by using an auxiliary laser rather than by a homogeneous heating of the whole part as problems of surface quality were encountered by Hagedorn et al. [START_REF] Hagedorn | Net shaped high performance oxide ceramic parts by selective laser melting[END_REF]. Different positions of the auxiliary laser will be firstly compared by cases #2 and #3. As the melt pool length with only principal laser is about 250 µm, ∆X 12 = 250 µm means the auxiliary laser is located at the melt pool tail. The auxiliary laser is activated only when it is between X = 0.3 mm and X = 2.55 mm in order to keep the same additional heating duration as illustrated by the red arrows in Figure

Single track without auxiliary laser

The mesh blockage strategy presented in the previous chapter provides a good representation of track shape. However, the mesh adaptation behind the melt pool leads to the coarsening of mesh elements. This effect results in the loss of precision in variable interpolation during mesh adaptation and it is more severe for element-wise stress and strain in the resolution of solid mechanics. Consequently, it is proposed to block the element when both the following conditions are satisfied:

• the element is not exposed to the laser radiation (the indicated powder at t = 2 ms in The stress distribution is shown in Figure 4.12 in x, y and z directions at t = 8 ms. The stress distribution appears to be in stable regime at this time. When the material under laser is heated, the decrease of density at higher temperature induces thermal expansion in this zone. This is inhibited by surrounding material, leading to local compressive stress. At the front of laser, material is in traction in direction y and z to inhibit this expansion. As a consequence, cracks may be formed by these tensile stresses, a part of them being later melted when the laser arrives. Note that this could justify the use of a second auxiliary laser working in front of the main one and acting as a preheater of the powder bed. Behind the melt pool, tensile stresses in x direction are generated during cooling as material has tendency to shrink but is inhibited by the surrounding. The stresses increase when temperature decreases as indicated along the formed track. Between 1000 and 2004 • C, stresses are below 1.5 GPa while they can atttain 3 GPa at T < 500 • C. In addition, σ xx is not homogeneous at different positions of remelted zone and it seems to be related to the track height. This will be discussed later by plotting the stress distribution profile. The stress distribution in y direction is more complicated. From the track surface to deep remelted zone, material experiences an alternation of zones in compression and traction. The stress magnitude of σ yy is much lower than σ xx . In z direction, the remelted zone is in traction-compression-traction. Particular attention should be paid to the high tensile stress at the beginning of track. The tensile stresses in x and y directions can explain the experimentally observed cracks in yz (Figure 1.8) and xz planes, respectively. Note that transverse cracks are more likely to occur in yz plane as σ xx is much higher than the tensile strength (240∼300 MPa). On the other hand,

STRESS DISTRIBUTION IN SLM

the delamination of remelted zone is related to tensile stresses in z direction. In Figure 4.13, the stress profiles are plotted along the X direction at initial powder/substrate interface in the median plane (horizontal white line in Figure 4.12). The local track height is also indicated with the associated scale on the right of the figure, together with a top view of the track (with powder). The maximum tensile stress in x and z directions is about 3.3 and 0.8 GPa, respectively while it is almost 0 in y direction along this profile. The increase of stress during cooling can be noticed. A remarkable finding is the correlation between local stress and the local track height (or track width). The stress peaks of σ xx , σ yy and σ zz correspond to the track necks (lower and narrower zone, so smaller section) and vice versa. It is interesting to see that peaks form at the same location independently of the stress component σ xx , σ yy and σ zz . This means that stress concentration can be generated in track necks. Consequently, track with uniform transverse section should be guaranteed from the point of view of geometrical precision and mechanical properties. 4.12) and track surface (with powder) in case #1. Note that the temperature increases along X and the laser is at X = 2.6 mm at this time t = 8 ms.

The stress distribution in depth is shown in Figure 4.14. The maximum tensile stress in x direction is located at the initial substrate/powder interface, similar for the maximum of σ zz . Compressive stresses of σ yy are found near the track surface, while the maximum tensile stress of σ yy is at the bottom of the remelted zone. This is in coherence with the result of Mukherjee et al. [START_REF] Mukherjee | An improved prediction of residual stresses and distortion in additive manufacturing[END_REF]. Hence, the cracks in xz plane may be initialized from the bottom of remelted zone (maximum σ yy ) to track surface. For the cracks in yz plane, it may be initialized at the initial position of substrate/powder interface (maximum σ xx ) and propagate to the track surface and bottom of remelted zone. It should be pointed out that the compressive stresses of σ yy near the track surface are not present in the results of Gusarov et al. [START_REF] Gusarov | Thermoelastic Residual Stresses and Deformations at Laser Treatment[END_REF] and Vastola et al. [START_REF] Vastola | Controlling of residual stress in additive manufacturing of Ti6Al4V by finite element modeling[END_REF]. This may be related to the track geometry, process parameters and material properties. The computation time is shown in Table 4.3. The total time is about 78.7 h for 4333 time steps until the extinction of the principal laser. The most time comsuming procedure is the mesh adaptation. Although it is done every 5 time steps, it represents still 24.9% of the entire computation time. The resolutions of fluid and solid mechanics are also time consuming, compared with the resolution of heat transfer, as 4 degrees of freedom per node are concerned in the resolution of fluid and solid mechanics. Further optimization of computing efficiency may be focused on the desynchronization of the resolution. For example, the solid 120 4.5. STRESS DISTRIBUTION IN SLM mechanics may be solved every specific number, N (N > 1), of standard time steps.

Single track with auxiliary laser

The use of an auxiliary coaxial laser (case #2) enlarges the heating zone around the principal laser, leading to a deeper and longer melt pool as shown in Figure 4.15. As the lateral zones close to the solidified track are heated, the cooling rate decreases and temperature iso-contours are extended as shown in Figure 4.16, especially in x direction. Due to the temperature gradient mechanism [START_REF] Mercelis | Residual stresses in selective laser sintering and selective laser melting[END_REF], the reduced thermal gradient results in lower tensile stresses in x direction at t = 8 ms compared with case #1 in Figure 4.12. The stress distribution at t = 12 ms after the extinction of all lasers shows clearly the development of high tensile stress during cooling. In addition, the maximum tensile stress is achieved at T < 500 • C. The auxiliary laser in case #3 is behind the principal laser by 250 µm. This configuration only affects the tail of the melt pool, rather than its front. Consequently, the melt pool shown in Figure 4.15 at t = 8 ms becomes longer (even longer than in case #2) but its depth is almost the same as in case #1. As the additional heating is imposed directly on the melt pool tail where solidification occurs, the thermal gradient and cooling rate are furtherly reduced compared with case #2. This can successfully decrease the stresses near the track surface. High stresses occur principally in the remelted zone below the initial powder/substrate interface. By contrast, the region after the extinction position of auxiliary laser remains relatively higher stresses as shown in Figure 4.17 Case #4 is assimilated to a CO 2 auxiliary laser with increased absorption while keeping the same heating position as in case #3. In this situation, the melt pool shape is similar to that in case #3 but with a slightly shorter tail (green contour in Figure 4.15). The higher absorption means that the energy input of auxiliary laser is limited to a smaller depth from the track surface. As a result, the temperature iso-contours in Figure 4.18 are shallower than those in Figure 4.17. Consequently, the tensile stresses in case #4 are higher than in case #3 comparing the red region in Figure 4.17 and 4.18 and a CO 2 laser is less efficient than the Yb:YAG laser to minimize σ xx built during cooling.

The stress distribution in depth is compared for cases #1-4 in Figure 4.19 along the profiles indicated by red vertical lines in Figure 4.12 and 4.16-4.18. All theses profiles are at the position where the surface temperature attains 500 • C. Cases #2-4 with auxiliary laser show a decrease of maximum tensile stress and case #3 has the lowest value. Hence, auxiliary laser with deep penetration appears to be more efficient to reduce the maximum tensile stress. Although the maximum σ xx is higher than the tensile strength in all cases, the use of an auxiliary laser deep penetration can reduce the crack risk more efficiently. In addition, it is better to impose it in the cooling zone at the tail of melt pool.

Conclusion

Solid mechanics in SLM process is modeled by an elasto-viscoplastic model. The powder bed is assimilated to a Newtonian fluid as we are not interested in its mechanical behavior. The advantage of this model is its generality and it can be easily applied to other materials. The FE implementation and numerical resolution are then detailed. The change of variable is presented in the context of mixed resolution of velocity/pressure. Its significant optimization on reduction of computation time is stressed.

Simulations are addressed to the effect of an auxiliary laser on the stress distribution. High tensile stresses in the scanning direction (x) can be generated in the melted zone behind the melt pool, leading to potential transverse cracks. The tensile stresses in transverse (y) and build (z) directions are responsible to potential cracks in longitudinal plane and delamination, respectively. They are much lower than in the scanning direction. In addition, high surface roughness can cause stress concentration, meaning more cracking risk. The use of an auxiliary laser can reduce the maximum stress. It is found that an auxiliary laser at the tail of melt pool is more efficient in decreasing the maximum stress due to the low thermal gradient in the cooling zone. The source type of auxiliary laser is also important as it is shown that auxiliary laser with deeper penetration (Yb:YAG) is better in the reduction of maximum stress than surface limited laser (CO 2 ). These results can be considered as guideline for the choice of auxiliary laser and its heating strategy in order to produce crack free parts.

Résumé

Dans ce chapitre, la comparaison entre simulations et expériences est présentée. Les expériences sont réalisées par dépôt d'un cordon dans différentes conditions de travail (P L , v L ). Ceci conduit à différents bains liquides et morphologies du cordon. Comme la source de chaleur est essentielle à la comparaison de la forme du bain liquide, les paramètres de la source sont calibrés à l'aide d'un modèle analytique, en minimisant la différence de dimensions du bain liquide entre la prédiction et l'expérience. Les dimensions du bain liquide simulé montrent une cohérence avec l'évolution expérimentale par rapport aux différentes conditions de travail (P L , v L ). D'autre part, la comparaison de la morphologie du cordon est plus difficile, principalement en raison de la présence d'évaporation, d'éclaboussures de liquide et de dénudation de poudre, qui ne sont pas considérées dans notre modèle. Cependant, pour les cas où E l est faible où ces phénomènes ne sont pas très prononcés, la simulation est représentative en termes de tendance de variation de hauteur du cordon. Un cas d'étude de l'effet de "balling" avec une vitesse de balayage plus élevée est également présenté. La formation de cordon interrompu sous la forme d'îlots est étudiée et son influence sur la distribution de la température est également discutée.

Le modèle est ensuite appliqué au dépôt multi-cordon avec différentes stratégies de balayage. Il est montré que la direction de balayage et la distance entres deux cordons successifs influencent la morphologie de la surface d'une couche. L'effet de la distance entre deux cordons successifs est le plus marqué, tandis que la direction de balayage affecte principalement la transition entre deux cordons successifs. L'influence de la direction de balayage sur la température et la contrainte est également discutée. Le balayage en zigzag permet d'obtenir une température plus élevée et d'avoir un recouvrement plus important du bain liquide dans la région de retour. La simulation montre également que le réchauffement peut diminuer les contraintes de traction dans les cordons précédents. La contrainte maximale qui se forme dans un cordon (à partir du deuxième) est également diminuée par un effet de préchauffage local par un environnement chaud (cordons précédemment chauffés). Ces prédictions numériques pourront être comparées à des expériences dans le futur, en particulier pour comparer l'influence des stratégies de balayage sur la morphologie de la couche et l'occurrence de la fissuration.

EXPERIMENTS

The modeling of heat transfer, melt pool dynamics and solid mechanics have been presented in previous chapters, together with sensitivity studies with respect to material properties and process parameters. This chapter aims at validating the numerical model by comparison with experimental results, in particular the shapes of melt pool and solidified track. First, the experimental setting is presented. The results of some tests are then used to calibrate the parameters of heat source model with the help of an analytical model. With these parameters, single track simulations are carried out to compare with experiments. In addition, the application of this model is extended to multi-track deposition to investigate the influence of scanning strategy on temperature evolution, layer surface morphology and stress evolution.

Experiments

The experimental part of the C ÉFAL É project is carried out at CdM in parallel to the present work, in the framework of the PhD thesis of Liliana Moniz Da Silva Sancho, supervised by Christophe Colin, Marie-Hélène Berger and Jean-Dominique Bartout. A brief introduction of this experimental work will be given in the following. More details can be found in [START_REF] Moniz | Laser Beam Melting of Alumina: Effect of Absorber Additions[END_REF]. The work consists of parametric studies for processing Al 2 O 3 regarding the laser power and scanning speed. These studies can help to establish a process window in which stable melt pools and track shapes can be obtained. These experimental results will be used for the validation of the developed model. Due to the transparent effect of Al 2 O 3 to Yb:YAG laser (see Section 2.3.1), it is difficult to measure the temperature distribution by thermographic camera as infrared radiation comes not only from the material surface, but also from the material depth. As a consequence, we will focus only on the comparison of melt pool and track shape for the moment.

Experimental configuration

The machine used is the PM100T from Phenix Systems. It uses a Yb:YAG laser (wavelength λ L = 1070 nm) with a power varying from 0 to 168 W. The dimension of its building chamber are diameter 100 mm and height 100 mm. A substrate is put at the bottom of the building chamber. It is prepared by cutting a relatively dense Al 2 O 3 cylinder (φ = 30 mm, 4% porosity) into a tablet with a thickness of about 5 mm. It should be mentioned that the existence of impurities in the substrate has significant influence on the energy absorption. Preliminary tests with laser treatment on the substrate show random fusion. As the absorption coefficient of pure Al 2 O 3 is very low, the Yb:YAG laser cannot melt it. However, the random presence of impurities with relatively higher absorption can increase the local absorption, leading to the fusion at these positions.

The Al 2 O 3 powder is mechanically crushed, leading to particle shape as shown in Figure 5.1, which is totally different compared with spherical powder particles obtained by atomization. The powder size distribution is measured with the laser diffraction method [START_REF] Slotwinski | Characterization of Metal Powders Used for Additive Manufacturing[END_REF]. It is described by the volume fraction function D V (50) = 14.9 µm, meaning that 50% (in volume) of powder particles are smaller than 14.9 µm. Processing Al 2 O 3 with Yb:YAG laser by SLM has additional difficulty compared with metals, due to its low energy absorption. Consequently, a certain amount (1wt% in the following cases) of carbon absorbers are added to increase the absorption. Figure 5.1 shows that the absorbers (average size of 42 nm) are heterogeneously distributed on the powder surface. Some particles are fully covered with carbon, while others are partially covered or not covered at all. The mixture of powder and absorber particles should be furtherly optimized. The powder is deposited on the substrate (Figure 1.2) with a layer thickness of 51 µm and porosity of 61%. Note that the layer thickness can be controlled by the displacement of the piston but the precise value is measured between the powder and the substrate surface after all tests have been realized as the substrate may be inclined in its initial state. On the other hand, the powder porosity may be different before and after layer deposition due to the slight tapping effect of the recoater during powder deposition. All tests are carried out within the same layer as indicated in Figure 5.2. For each given laser power of 42, 84, 126 and 168 W, 5 tests with scanning speed of 50, 100, 200, 400 and 800 mm • s -1 have been done one beside one. In the following, they will be noted like P42V50 for the case with P L = 42 W and v L = 50 mm • s -1 , and similar for others. The powder is then removed, leaving the whole sample with 20 solidified tracks attached to the substrate. The track surface is observed with a Scanning Electron Microscope (SEM) and measured by an optical profilometer. The sample is also cut into transverse section to measure the melt pool dimensions.

Experimental results

Effect of scanning speed on denudation

The colored surface height of powder and solidified track after single track scanning without any post-treatment is shown in Figure 5.3 for the cases with P L = 84 W. The height increases from green color to yellow. One can notice that there is a low region beside the tracks. At low scanning speed (e.g. 50 mm • s -1 ), this region may be explained by the shrinkage in width under the effect of surface tension. However, when the scanning speed increases, this region is significantly enlarged while the track becomes narrower. Hence, there must be powder denudation in these cases. In cases P168V50 and P168V100, the linear energy is so high that the track is damaged because of the internal constraints. This may be caused by the keyhole effect discussed in Section 1.3.2. Under this effect, liquid spatters may be formed due to the evaporation pressure, leading to a destroyed track. In addition, track irregularity with width variation can be remarked in several cases, as in P126V100. This is probably due to the stochastic nature of the powder along the lateral sides of the track. This has eventual consequence on the neighboring tracks as the wetting of melt pool is different. In all cases, there are periodical cracks in the yz plane caused by high tensile stress σ xx . This has been demonstrated in our previous simulations in Section 4.5. Numerical results show that tensile stresses can be decreased by using an auxiliary laser with appropriate heating strategy. Hence, cracks can be hopefully eliminated, being a topic to be investigated in the future.

EXPERIMENTS

Melt pool shape

The sample is cut in transverse section to observe the shape of melted zone as shown in Figure 5.5. Due to the shrinkage of alumina during its solidification, the melted zone is stripped from the non-melted substrate, leading to a crack boundary corresponding to the melt pool. Similar to the track width, the melt pool becomes deeper (note the different scales in Figure 5.5) with the increase of laser power and the decrease of scanning speed. For some cases with high linear energy (e.g. P168V50 and P126V50), the melt pool appears to be 100 µm µm is used in all other cases. Melt pool shape of cases in red rectangles with linear energy E l between 0.2 and 0.9 J • mm -1 is considered more suitable for the construction. Note that the melt pool shape for a given condition may vary with different cutting positions.

like a bolt, with a very deep penetration into the substrate. By contrast, at a high scanning speed of 800 mm • s -1 , the substrate is almost not melted. For the process, linear energy E l between 0.2 and 0.9 J • mm -1 is considered as the suitable process conditions for desired melt pool shape, as indicated in red rectangles. In addition, the crack boundary is much more visible in cases with higher linear energy E l . This is caused by the higher tensile stresses on the boundary with higher E l , which can easily provoke the cracks. Note that the melt pool shape may vary with different cutting positions. In fact, residual porosities are present in most cases but may be not seen in the chosen cutting section in Figure 5.5. These porosities are probably due to the rapid solidification of liquid with entrapped gas. Melt pool dimensions, including the width, depth and height, are defined in Figure 5.6. The width and depth are measured and plotted in Figure 5.7 with a systematic measurement uncertainty of 5 µm. The increase of with and depth with the increase of laser powder and the decrease of scanning speed are systematic. However, the cases with P L = 168 W and v L ≥ 200 mm • s -1 show no increase of melt pool depth, compared to cases with P L = 126 W.

MEASUREMENT AND CALIBRATION OF PARAMETERS DEFINING THE HEAT SOURCE

For all cases with v L = 800 mm • s -1 , the penetration of the melt pool into the substrate is very small. There may be even no remelting of substrate in the case P42V800.

Measurement and calibration of parameters defining the heat source

The heat source model in Eq.2.22 depends on three parameters: the reflection coefficient R, the effective interaction radius r int and the local absorption coefficient α. Among them, the local absorption may vary with respect to the temperature and the material state, making the calibration of the heat source model more difficult. In the following, the measurement of R by an optical method will be firstly presented. However, the parameters r int and α are difficult to measure directly. Therefore, they are calibrated by an inverse approach using numerical simulations and an analytical model detailed hereafter.

Reflection

Laser entrance The reflection coefficient at the surface of the powder bed is measured by an optical method using an integrating sphere as shown in Figure 5.8(a). This instrument consists of a hollow spherical cavity with small holes for entrance of laser source and exit ports. The sample is put at the position facing the laser entrance and a side-detector is used to measure the light intensity. A gain unit is cabled with the detector for amplification and gives the final output signal (in voltage). The most important part is the inner surface coated with diffusive and reflective material, such as barium sulfate, diffused gold, depending on the intended application and light wavelength. This material provides uniform scattering effect as illustrated at point A in Figure 5.8(b), leading to equally distributed light intensity in all directions.

The measurement of reflection begins with a reference sample whose reflection coefficient R ref is known a priori. With a given laser power, the output signal can be noted as S ref .

Then the same measurement is done with the same laser power for the sample to be analyzed. Assuming that the scattering on sample surface is also uniform in all directions, the reflection coefficient R of the sample can be simply calculated by:

R = S S ref R ref (5.1)
where S is the output signal of sample. Note that possible errors come principally from the hypothesis of homogeneous scattering of the sample surface. The measurement of reflection coefficient R is carried out with different laser power and scanning speed. The result is plotted in Figure 5.9 as a function of linear energy E l . For some conditions, such as E l = 0.105, 0.21 and 1.68 J • mm -1 , two tests are conducted. The error bar is deduced from the temporal evolution of measured signal. Note that what is measured is the average reflection between the powder and melt pool as the laser irradiates both of them. It shows that the reflection decreases with the increase of linear energy. As indicated, low E l (top left) means that the laser irradiates mostly the powder. By contrast, with high E l (top right), most of the irradiation is imposed on the melt pool. Consequently, the results clearly show that powder has higher reflection (about 8%) than liquid (about 2%). In the following, we will not distinguish the reflection of powder and liquid in order to use the analytical model CHAPTER 5. EXPERIMENTAL VALIDATION AND OTHER APPLICATIONS frame whose origin is attached to the start point of the laser:

φ 0 (x, y, t) = 2P L πr 2 int exp -2 (x -v L t) 2 + y 2 r 2 int (5.4)
The main idea of this analytical model is to calculate the energy received by a local material point from laser. The energy increase ∆h v of points at melt pool boundary Z l should correspond to the fusion enthalpy ∆h v,f :

∆h v (x, y, Z l , t) = ∆h v,f (5.5)
This can be distinguished into two regimes: (1) Z l = 0, which is before the initialization of fusion on surface; (2) Z l > 0, where the melt pool is formed and in development.

Regime 1: Z l = 0

This case corresponds to the minimum condition of fusion. With the assumption of no heat exchange and considering that the laser scans from t 1 = 0 to t 2 = t, ∆h v (x, y, 0, t) of a point located at the mid-trajectory of laser (x = v L t/2) can be obtained by integrating the heat source (second case of Eq.5.3) from 0 to time t. This results in (see more details in Appendix):

∆h v (x, y, 0, t) = α s (1 -R)E l √ 2πr int exp -2 y 2 r 2 int erf √ 2 v L t -x r int -erf √ 2 0 -x r int (5.6)
where erf (x) is the error function defined by erf (x) = 2 √ π

x 0 e -u 2 du. Consequently, the fusion condition is:

∆h v,f ≤ ∆h v (x, y, 0, t → ∞) = 2 π α s (1 -R)E l r int exp -2 y 2 r 2 int (5.7)
The fusion initialization time t f (x, y) can be deduced by ∆h v = ∆h v,f in Eq.5.6.

Regime 2: Z l > 0

Once the fusion condition in Eq.5.6 is satisfied, the melt pool can form and develop. As there is no heat exchange, the melt pool will attain its maximum size Z f,max when t → ∞. In the following, we aim to deduce the expression of Z f,max .

The derivation of Eq.5.5 with respect to time leads to:

∂∆h v ∂ z (x, y, Z l , t) ∂Z l ∂ t (x, y, t) + ∂∆h v ∂ t (x, y, Z l , t) = 0 (5.8)
Taking Z f,max as unknown, the melt pool shape can be deduced by solving this differential equation. Due to the assumption of no heat exchange, the received energy of local point from

MEASUREMENT AND CALIBRATION OF PARAMETERS DEFINING THE HEAT SOURCE

laser is conserved, meaning:

∂∆h v ∂ t (x, y, Z l , t) = qL (x, y, Z l , t) (5.9)
On the other hand, Defillon et al. [START_REF] Defillon | Modèle du bain de fusion de SLM sur céramique avec disparition de l'agent absorbant restreinte au bain[END_REF] showed that:

∂∆h v ∂ z (x, y, Z l , t) = -α s ∆h v,f (5.10)
Inserting Eq.5.3, 5.9 and 5.10 into Eq.5.8 leads to:

∂Z l ∂ t = 1 -R ∆h v,f φ 0 exp (-α l Z l ) (5.11) or ∂ ∂ t exp (α l Z l ) α l = 1 -R ∆h v,f φ 0 (5.12)
Note that both Z l and φ 0 are functions of x, y and t. By integrating Eq.5.12 from the fusion initialization time t f to infinity, Z f,max can be finally obtained:

Z f,max (y) = 1 α l ln   1 + α l   2 π (1 -R)E l r int ∆h v,f exp -2 y 2 r 2 int - 1 α s      (5.13)
This expression gives a full description of the melt pool boundary in the transverse plane. Hence, the melt pool depth H zr and width w app can be easily obtained:

H zr = Z f,max (0) = 1 α l ln   1 + α l 2 π (1 -R)E l r int ∆h v,f - 1 α s   (5.14) Z f,max w app 2 = 0 → w app = 2r int 1 2 ln α s (1 -R)E l r int ∆h v,f - 1 4 ln π 2 (5.15)
The melt pool width and depth predicted by this analytical model can be considered as first approximations. It can be noticed that they depend on the laser scanning speed and power only through the linear energy E l . In real process or simulation, the melt pool dimensions will be different due to the presence of diffusion and convection.

Effective interaction radius and absorption

The effective interaction radius and absorption coefficient have direct consequence on the melt pool shape. Hence, the melt pool dimensions like width and depth can be used to calibrate r int and α. They may be obtained by minimizing the width and depth differences between experiments and calculations under different process conditions. Proceeding with calculated values obtained by finite element simulations would require a group of simulations for each given couple of (r int , α s , α l ), which is very time consuming. Instead, it is proposed to firstly estimate these parameters by minimizing the errors between experiments (superscript exp) and calculated values obtained by previous analytical model (superscript ana) as follows:

(r int , α s , α l ) = arg min (r int ,αs,α l ) (P i L ,v i L )∈G   H i,ana zr H i,exp zr -1 2 + w i,ana app w i,exp app -1 2   (5.16)
where G is the set of process condition couples (P L , v L ) used for the parameter calibration. Note that each of the two contributions of error function is normalized, in order to give the same weight to depth and width. Eq.5.16 is solved numerically through the software Mathematica.

As discussed in Section 5.1.2, high or low linear energy leads to keyhole or absence of melt pool, respectively. Consequently, only process conditions (0.2 < E l < 0.9 J • mm -1 , cases in red rectangles in Figure 5.5) with suitable melt pool shape are used to calibrate the heat source parameters. In addition, all tests with P L = 168 W are not used as the melt pool depth is almost the same to P L = 126 W with same v L , which is not normal. Hence, only the process conditions in Table 5.1 are used for the calibration, which consists the input set G in Eq.5.16. The simulated system is similar to that in Section 3.6, but with a layer thickness of 51 µm and porosity of 61% of powder. The first calibration with Eq.5.16 gives r int = 73.1 µm, α s = 8.13 mm -1 and α l = 6.1 mm -1 , which results in predicted melt pool dimensions as shown by the dashed lines with cross markers in Figure 5.11. Using now these calibrated heat source parameters, the melt pool dimensions obtained by FE simulations are also given by solid lines with circle markers in Figure 5.11. Note that balling effect is observed in simulation for the case P42V200, so the melt pool dimensions are not plotted.

For the cases with same linear energy, the analytical depths (or width) are the same according to Eq.5.14 (or Eq.5.15). Although the numerical FE simulations are conducted with parameters calibrated with the help of the analytical model, there is remarkable difference between analytical prediction and simulation. This difference comes principally from the and supposing that this difference is the same under a fixed process condition while change slightly the heat source parameters, we can estimate a set of calibrated parameters that would be more appropriate to the use of the FE numerical simulation. Hence, it is proposed to calibrate heat source parameters by solving:

(r int , α s , α l ) = arg min (r int ,αs,α l ) (P i L ,v i L )∈G   H i,ana zr -∆H i zr H i,exp zr -1 2 + w i,ana app -∆w i app w i,exp app -1 2  
(5.18) Note that ∆H i zr and ∆w i app are evaluated with the first analytical prediction and simulation. This second calibration with Eq.5.18 leads to r int = 64.6 µm, α s = 8.98 mm -1 and α l = 3.88 mm -1 , thus mainly showing an adaptation of the absorption in the liquid. These values will be used in the following comparison between FE simulations and experiments.

Comparison of melt pool and track shape

In this section, simulations are carried out with only heat transfer and fluid mechanics as the objective is to compare the melt pool shape. The configuration of simulation is similar to those in Chapter 3, but with a powder layer of 51 µm and powder porosity of 61%. In addition, the reflection and calibrated parameters are used in the following cases.

Melt pool shape

The experimental results are overlaid with the melt pool shapes obtained by simulations as shown in Figure 5.12. Note that the melt pool shape may vary with different cutting positions. The simulated melt pools have the same shape as experiments for low E l . However, for the cases with deep penetration, the evolution of section shape in depth is different as indicated in the blue rectangle for the case P42V50. .12: Comparison of melt pool shape for different process conditions. The experimental results are overlaid with simulated melt pool contours (red profiles) with the same scale. Note that the melt pool shape for a given condition may vary with different cutting positions.

COMPARISON OF MELT POOL AND TRACK SHAPE

Quantitative comparison of melt pool width and depth is shown in Figure 5.13 for different laser power and scanning speed. Compared with the first simulation results in Figure 5.11, heat source parameters deduced with correction in analytical model result now in less differences between experiments and simulations. Among them, the cases P126V200, P126V400, P84V200 and P84V400 show good coherence with experiments, in both width and depth. The principal differences are found for the cases with P L = 42 W. The melt pool is smaller in simulation than in experiment both in width and depth. The procedure of parameter calibration may be done again to furtherly minimize the difference.

Track shape

The track height measured by an optical profilometer is compared with simulations in Figure 5.14 for the cases P42V200 and P84V400. Experiments show clearly the height variation in the scanning direction. This is also predicted by the simulation as a result of melt pool dynamics. For the case P42V200, the track height map of simulation shows coherence with the experiment. The track width variation can be also noticed by the regions with attachment of the lateral sides of track to the partially melted powder as indicated. The track is also higher in these regions, meaning larger cross sections at the same time. Between some of these regions, necking occurs, which is the potential position of track interruption. The comparison of track height profile is then taken along the median plane. Only the track between X =0.5 mm and 2.5 mm is plotted for the simulation as this is the stable regime. High frequency variation can be noticed in experiment with sharp peaks or valleys. In some positions, the track height is even negative, meaning the absence of track. It is probably caused by liquid spattering, which is not integrated in the model. If we neglect this effect and the high frequency variation, simulation matches well for the experiment both in the average height and variation amplitude. For the case P84V400, the track height map and profile indicate more regions lower than 20 µm, so more spattering exists in this case. Apart from these regions, the simulation shows the tendency of height variation although the average height in simulation is higher than in experiment.

Track regularity

In SLM, the balling effect is a severe situation of track irregularity, as introduced in Section 1.3.2 (Figure 1.7). It results in poor surface quality and mechanical performance as the deposited layer is not well attached to the previous consolidated layer. This happens when the ratio of melt pool length over cross section attains the Rayleigh-Plateau limit [START_REF] Gusarov | Heat transfer modelling and stability analysis of selective laser melting[END_REF]. This condition can be met with high scanning speed. Hence, additional case is carried out to study 5.3. COMPARISON OF MELT POOL AND TRACK SHAPE this effect. Figure 5.15 shows the appearance of balling with P L = 84 W and v L = 600 mm • s -1 , which is higher than the scanning speed of previous cases with P L = 84 W. The track is not continuous any more, but broken into isolated islands, just like the balling effect in Figure 1.7 at v L = 500 mm • s -1 . These islands seem to be periodical along the scanning trajectory, except the fifth, which is connected to the next. The height map and profile show that the maximum height of islands is about 50 µm, which is close to the layer thickness. The zone between isolated islands is at the same level of the initial substrate surface. This means that the liquid above the substrate is cleaned up and gathered together by surface tension to form these islands. In order to understand the formation of such isolated track islands, a series of snapshots is shown in Figure 5.16 with black contours corresponding to the melt pool. At t = 1 ms, the principal melt pool (on the right-hand side) is developing by the fusion of powder. When it becomes longer at t = 1.1 ms, it tends to shrink under the effect of surface tension at t = 1.11 ms, leading to a neck in the contact region between liquid and powder. As the energy is blocked in this small melt pool and heat extraction by substrate is limited by the small contact area, the maximum temperature at t = 1.1 ms rises up to 2896 • C, even higher than the value (∼ 2370 • C) in case P84V400 with higher linear energy. Note that this is close to the boiling point (2970 • C), meaning more risk of evaporation and liquid spattering. The neck is then broken at t = 1.13 ms. The melt pool retracts in -X direction at t = 1.2 ms, forming the third melt island. The fusion of powder continues to develop the principal melt pool. This phenomenon is periodical and a series of melt islands is left behind the laser. These melt islands solidify and result in the final track in the form of track islands. An exception occurs for the principal melt pool at t = 1.8 ms. It is not totally broken and keeps connection with the next one (t = 2 and 2.3 ms), leaving two solidified islands connected (t = 2.3 ms). Behind the laser, the cooling is uneven as temperature is higher in islands than in the region between them as shown at t = 4.4 ms. Consequently, temperature peaks appear along the profile just below the initial substrate surface, leading to alternate sign of temperature gradient . This can eventually affect the stress distribution. Comparing the literature result in processing 316L [START_REF] Li | Balling behavior of stainless steel and nickel powder during selective laser melting process[END_REF], numerical simulation is able to reveal this phenomenon, although the applied materials are different.

Application to multi-track deposition

All previous simulation cases are based on single track deposition. However, the developed model is applicable to the simulation of multi-track deposition, which is representative for the layer deposition. In this section, we aim at studying the influence of scanning strategy on thermal evolution, layer surface morphology and stress evolution during multi-track deposition.

Simulation configuration

Previous case P84V400 shows that a height of 0.949 mm (leading to total height of material equal to 1 mm with powder layer thickness 51 µm) for the substrate may be not necessary as the melt pool depth is limited to about 50 µm and the bottom is not highly thermally affected. Consequently, in the following multi-track simulations, a substrate with height of 0.449 mm is used, with a layer of powder of 51 µm deposited on it, as shown in Figure 5.17. Above them, there is still the gas domain with a height of 0.1 mm. Different scanning strategies in multi-track deposition can be used by the combination of scanning direction and hatch distance ∆y (see Figure 5.18). These two factors are investigated in the following by 3 cases listed in Table 5.2 with P L = 84 W, v L = 400 mm • s -1 and the same heat source parameters used in the previous single track validation. For each case, 3 passes (numbered from #1-3) are carried out as illustrated in Figure 5.18 for unidirectional scanning. For the zig-zag scanning, pass #2 is in inverse direction from right to left. ∆y = 50 µm for cases #I and II is chosen to be close to the track width (Figure 5.13). The transition between the end of previous pass and the beginning of the current one is considered to be instantaneous. Finally, case #III is the same as case #I but reducing the hatch distance ∆y to 40 µm.

Simulations are conducted with the full algorithm presented in artificial viscosity of liquid, gas and powder to 10 5 Pa • s in the following simulations. The number of elements evolves from about 1 320 000 to 2 080 000 from the beginning to the end as the surface of all three tracks are refined.

Layer surface morphology

The evolution of the layer surface for each case is illustrated in Figure 5.19 at the end of each pass, by the colormap of layer height with respect to the initial substrate surface. Pass #1 shows a result close to the single track case P84V400 (Figure 5.14(b)) with slightly changed track surface height. Differences between different scanning strategies can be noticed from pass #2. For case #I, there is a clear boundary in the shape of groove, between passes #1 and #2. In case #II with zig-zag strategy, two regimes can be distinguished. The first regime corresponds to the transition region (indicated in the red rectangle) when the laser turns back from pass #1 to #2. Two tracks are then connected without visible boundary. The boundary appears again out of this region and the surface looks similar to case #I. The decrease of the hatch distance to 40 µm in case #III consists in increasing the overlap between the successive tracks. This leads to a boundary almost invisible: the two first tracks appear like a wider single track. Note that the highest regions in cases #I and II are found at the same X positions. This means that the morphology of the previous solidified track may have an influence on the morphology of the current one. In addition, for all cases, one can see that the non-overlapped lateral side of pass #2 is lower than that of pass #1 Figure 5.19(b). However, the two lateral sides of track #1 are almost symmetric (Figure 5.19(a)).

APPLICATION TO MULTI-TRACK DEPOSITION

During pass #2, the melt pool seems to be pushed towards the 1st solidified track, leading to this difference. For track #3, there is still a boundary with track #2 in case #I and even more visible, while the boundary occurs partially again in case #II. The difference is in case #III, where the boundary appears, rather than being smoothly connected. All three cases show the decrease of surface height fluctuation from track #1 to #3, which again approves the influence of previous solidified tracks on the surface morphology. The profile evolution of the gas/material interface (ψ = 0) in the cross section X = 1.5 mm is plotted in Figure 5.20, taken at the end of each pass. Along the lateral sides of the developed track, there is still some powder in contact with the track as indicated in Figure 5.20(b). Due to weak overlap in cases #I and #II and the concave shape of track section, the groove between tracks #1 and #2 can be clearly observed, while it does not appear in case #III with ∆y = 40 µm, leading to a smoother surface. The height of track #2 is lower than #1 in cases #I and #II, but they are almost the same in case #III. Track #3 is more expanded in the width direction y and the groove appears in all cases. In fact, under the surface tension effect, the melt pool tends to be attracted by the previously solidified track. The slight nonphysical displacement of interface is noticed as shown in Figure 5.20(c). It has tendency to decrease the interface curvature. This is caused by the mesh adaptation. Although elements are blocked when the track is solidified, they are unblocked when they are exposed to the laser irradiation (numerically truncated to 1.5r int from the laser axis). In order to better preserve the solidified track surface, this blockage strategy may be optimized. The variation of layer surface results in the problem of mass conservation. The evolutions of mass deviation (see definition in Eq.3.55) are plotted for all three cases in Figure 5.21. The mass of material domain increases and the deviation tends to a stable value between 6% and 8%.

Figure 5.21: Mass deviation in the three simulation cases of multi-track deposition. Note that pass #1 is the same for all, so the 3 curves are overlaid between t = 0 and 6.5 ms.

Temperature and stress distribution

The layer surface morphology is a direct result of the evolution of melt pool shape. Temperature distribution and melt pool shape are thus compared in Figure 5.22 for cases #I-III, taken at several special instants. At the end of pass #1 (t = 6.5 ms, same for all cases), the beginning of track #1 has been cooled to around 200 • C. The adjacent powder is heated only in a limited width due to its low conductivity. At the beginning of pass #2 at t = 6.9 ms, the laser affects the already cooled beginning of track #1 in case #I, while the end of track #1 is still solidifying, leading to separated melt pools. On the other hand, with the zig-zag scanning in case #II, the laser turns back at the end. As a consequence, one big continuously evolving melt pool is obtained. This explains the formation of two connected solidified tracks in the turning back region mentioned in the previous section (red rectangle in Figure 5.19(b)). In this case, material is irradiated by the laser for a longer duration. Hence, the maximum temperature quickly arises to almost the boiling point (2970 • C), meaning potential risk of evaporation. Despite of the difference between cases #I and #II, when the laser arrives at the middle of pass #2 at t = 9.75 ms, both cases #I and #II appear with a same melt pool shape and they both have the stationary maximum temperature near to 2440 • C, as for pass #1. In addition, they have almost the same overlap with track #1, which explains the same groove position between tracks #1 and #2 for case #I and II in Figure 5.20(a) and (b), respectively. Comparing the melt pool shape at this time with the one at t = 6.5 ms, it has a relatively sharper tail. This is due to the heat extraction by the contacted track #1, which has a higher conductivity than powder. In case #III with reduced hatch distance ∆y = 40 µm, the melt pool is closer to the first scanning line (white line) than in case #I. It can almost wet the top of track surface, leaving the two solidified tracks appearing as in Figure 5.20(c). The beginning of pass #2 in case #III is similar to that in case #I and pass #3 is similar to pass #2, so they are not presented. In multi-track deposition, the material undergoes several cycles of heating and cooling. This thermal history leads to different stress distribution than in single track deposition. Figure 5.23 shows the axial stress distribution σ xx in the XZ plane at the position Y = 0.25 mm (scanning line of pass #1). Stresses develop during solidification and high tensile stresses near to 3.8 GPa are reached in the first half of track #1 at t = 6.5 ms (the end of pass #1), which is much higher than the tensile strength (240∼300 MPa). During pass #2 CHAPTER 5. EXPERIMENTAL VALIDATION AND OTHER APPLICATIONS (t = 6.5 ∼ 13 ms) and for both cases #I and #II, the previously developed tensile stresses are progressively eliminated due to the reheating (by the laser passed just beside) as shown at t = 9 ms. In addition, as the surrounding is hot (meaning low thermal gradient), the stress development is limited.

Temperature history

Heating and cooling cycles can be investigated by the temperature history at local points as shown in Figure 5.24 for the case #I. During pass #1, points 1-3 are successively and quickly heated near to 2190 • C, just above the liquidus (T l = 2104 • C). The heating rate is estimated to be 5.3 × 10 6 K • s -1 (averaged between the beginning and the time at which the maximum temperature is reached). They are then relatively "slowly" cooled, compared with the high heating rate. However, the mean cooling rate is about 3. From this temperature, it is heated and melted at the beginning of pass #3, with a maximum temperature slightly higher than the value attained by points 1-6. In this unidirectional scanning, all points have the same period of heating and cooling cycle.

In addition, the cooling rate at points 3 and 6 are higher than at other points as the laser switches instantaneously to another sides after heating these two points. Case #II in Figure 5.25 with zig-zag scanning shows different temperature evolution at points 1, 3, 4, 6, 7 and 9, who are located in the turning back region. They undergo a longer duration of heating and cooling. In fact, the additional cooling does not largely decrease the temperature as at the point 1 during pass #2. However, additional heating largely increases the maximum temperature up to 2500 • C at points 6 and 9. Points 3 and 4 are reheated just after a short period of cooling as they are again irradiated when the laser begins passes #2 and #3, respectively. For the other points (2, 5 and 8) at the middle of each pass, they go through almost the same temperature history as in case #I. The comparison between Figure 5.24 and 5.25 shows that different temperature histories are located in the turning back region. Hence, for multi-track deposition with long pass, there may be no significant difference of temperature history in the center.

Stress history

The stress evolution at each point can be studied in addition to the temperature evolution. As the axial stress component σ xx is the most detrimental to the cracks (in Y Z plane) observed in experiments (see Figure 5.4), it is firstly investigated in details as shown in

APPLICATION TO MULTI-TRACK DEPOSITION

for case #I. At the beginning of pass #1, point 1 is in contraction in x direction due to the thermal expansion as the laser is just behind this point. It is interesting to see that at this time, point 4 is also in contraction in x direction while point 7 is in tension, meaning a change of sign of σ xx along Y direction in the front of the laser. σ xx becomes to be 0 when point 1 is later melted (see temperature evolution in Figure 5.24). During the cooling, the tensile stress at point 1 increases rapidly and then stay in a stable regime around 3 GPa. Point 2 attains a maximum tensile stress higher than points 1 and 3, with a value near to 3.6 GPa. For all points 1-3, reheating when the laser passes beside them decreases quickly the stress. After that, the stresses increase much more slowly, mainly due to the lower thermal gradient. At points 4-7 (not representative for points 8 and 9), σ xx increases also slowly as their surrounding is already heated by pass #1. In addition, their maximum tensile stress is below 3 GPa due to the same origin in thermal evolution. Note that this also proves the effect of preheating on the possible elimination of cracks through the smoothing of higher stress levels. A strange evolution is observed at point 3 with oscillations appearing with regular envelop. In fact, it is due to the mesh adaptation as the transition between plateau and valley (see the zoomed rectangle) is always at the moment of mesh adaptation. This is not observed in temperature evolution. This difference comes from the fact that the temperature is a nodal variable (P1) while the stress is based on element (P0). A point located at the boundary of two neighboring elements may be easily alternated between the values of the two elements (P0 variable) due to the numerical precision. By contrast, a nodal variable takes always linear interpolation. Although cracks in XZ plane are not observed in experiments (Figure 5.4), they can be provoked by the stress component σ yy . Hence, the local evolution of σ yy at points 1-3 are also plotted in Figure 5.27(a) for cases #I. For pass #1 at about t = 3 ms, the laser arrives just behind the point 2. Material at point 2 is in compression in direction x (σ xx ≈ -1 GPa in Figure 5.26) while it is in tension in direction y with σ yy ≈ 1 GPa. This value is higher than the tensile strength, so cracks may occur. However, they are then eliminated due to the fusion of material. During pass #2, the stress state is alternated from tension to compression at point 1. Tension stress is developed at point 2 and 3 until the laser passes beside them. During pass #3, stress at point 1 is again alternated, from compression to tension, with a stable σ yy about 500 MPa. It is higher than the tensile strength, but cracks in XZ plane are not observed in experiments. σ yy at point 2 becomes negative at the beginning and then increases near to zero. At point 3, σ yy increases after the laser passes beside it. One can notice that the stress histories at points 1-3 are totally different, mainly due to their different temperature histories and the temperature state of the surrounding.

APPLICATION TO MULTI-TRACK DEPOSITION

The influence of the scanning direction on the stress component σ xx is studied in Figure 5.28 for case #II to compare with Figure 5.26. It is similar to its effect on temperature as shown when using a zig-zag scanning. Points 2, 5 and 8 have the same evolution as in case #I and differences are always restricted to points in the turning back region. At points 1, 3 and 6, material endures a longer duration at high stress than at other points. The evolution of σ yy with different scanning directions can be also compared in

Conclusion

In this chapter, the comparison between simulations and experiments is presented. Experiments are carried out by single track deposition under different process conditions (P L , v L ). This leads to different melt pools and track morphologies. As the heat source is essential to the comparison of melt pool shape, its parameters are calibrated through an analytical model, by minimizing the difference of melt pool dimensions between prediction and experiment. The simulated melt pool dimensions show an agreement with experimental evolution with respect to different process conditions (P L , v L ). The melt pool width and depth match well for cases with P L = 84 and 126 W, while differences exist for P L = 42 W. On the other hand, the comparison of track morphology is more difficult, principally due to the presence of evaporation, liquid spattering and powder denudation, which is not taken into account in our model. However, for cases with low E l in which these phenomena are not very pronounced, the simulation can still show similar track height variation. A study case of balling effect with higher scanning speed is also presented. The formation of interrupted track in the form of islands is investigated and its influence on temperature distribution is also discussed.

The model is furtherly applied to multi-track deposition with different scanning strategies. It is shown that both the scanning direction and hatch distance influence the layer morphology. Among them, the effect of hatch distance is more important, while the scanning direction affects mainly the transition stage between two successive tracks. The influence of scanning direction on temperature and stress are also discussed. Zig-zag scanning results in higher temperature and more melt pool overlap in the turning back region. Simulation also shows that reheating can decrease tensile stresses in previous tracks. The maximum stress forming in a track (track #2 and followers) is also decreased due to a local preheating effect by a hot surrounding (previously heated tracks). This application may be compared with experiments in the future, especially to compare the influence of scanning strategies on layer morphology and cracking occurrence.

Conclusions and perspectives

Conclusions

During this research activity, a numerical model dedicated to the simulation of SLM process applied to ceramics has been proposed in the context of the project C ÉFAL É . Although the scope of the project C ÉFAL É regarding the selection of ceramics materials was broader than the sole alumina, this study offered an investigation of SLM process from a numerical point of view. In the following, we summarize the main achievements of our work according to different aspects.

• Modeling method

This numerical model is developed at the track scale using FE method. Continuous powder bed is assumed without the representation of powder particles. The choice of this scale is based on the consideration of reasonable computation charge and the physical phenomena to be modeled, mainly heat transfer, melt pool dynamics and stress evolution during track deposition. The level set method is used to model a two-domain system and track the evolution of the gas/material interface induced by the condensation from powder to dense matter and the hydrodynamics in the melt pool. A full description of the multiphasic material domain is used, offering the details of phase evolution during heating and cooling. The model takes advantage of anisotropic mesh adaptation available in Cemef.

• Heat transfer

In the modeling of heat transfer, a heat source model based on Beer-Lambert law has been implemented, considering the local absorption of material. The energy conservation equation is solved with the temperature-based solver previously developed in Cemef. It takes into account the thermodynamics of material and can be coupled with solidification path. The non-linear relation between temperature and enthalpy is solved by Newton-Raphson method. It offers full access to the material state, including temperature, enthalpy, phase fraction, density, etc.

The interaction between laser and material is the driving force of all physical phenomena. Hence, the heat source model describing the interaction is very important. We proceeded to analytical calculations (Figure 2.10) showing that it is very sensitive to the effective interaction radius r int and material absorption α. The latter may vary with different materials and simulation results (tests in Section 2.6.3) have shown its signif-CHAPTER 5. CONCLUSIONS AND PERSPECTIVES icant influence on the temperature distribution and melt pool shape. For the tested process parameters (tests in Section 2.6.4), the increase of laser spot size or scanning speed leads to shorter, narrower and shallower melt pool (or temperature iso-contours), while the laser power has contrast effects. It should be noted that with the low absorption of ceramics, simulation results show deep penetration of heat source, rather than being limited to the material surface as for metals.

• Melt pool dynamics For the modeling of melt pool dynamics, the important surface tension and Marangoni forces are integrated. The surface tension is implemented with a semi-implicit formulation, providing a more stable resolution. The mesh is crucial in our model and the mesh adaptation is detailed. It includes the metric construction based on error estimation, the intersection of metrics and the blockage of elements. Compared the fixed grid usually used in many literature results, the constructed anisotropic mesh in our model reduces efficiently the computation charge.

Fluid dynamics is important for the melt pool stability and track regularity. In addition, it affects the temperature distribution via convection. The formation of liquid droplets is observed in simulation as the melting powder tends to spheroidize under surface tension. The collapse of these droplets causes the surface oscillation of melt pool, resulting track surface roughness after solidification and possible balling effect under surface tension. Another important contribution comes from the Marangoni effect. It creates convection circles depending on the sign of ∂γ/∂T (Figure 3.22). This convection redistributes the heat and changes the shape of melt pool. One interesting finding is the smoothing of track surface with negative ∂γ/∂T and opposite effect with positive ∂γ/∂T . However, the rather high viscosity of ceramic materials (higher than metal) tends to mitigate the melt pool dynamics. The combined effect of surface tension and scanning speed is also an interesting point. High scanning speed means high solidification speed. Consequently, melt pool oscillations can be quickly "fixed" by solidification. Our simulations show high track irregularity (Figure 3.22-3.24) and even the occurrence of balling effect(Figure 5. 15-5.16).

• Solid mechanics

The modeling of solid mechanics is based on a previously developed elasto-viscoplastic solver for welding simulation. In the context of LS modeling for powder-based SLM process, we are only interested in the mechanical behavior in solid. All the liquid, gas and powder are considered as incompressible Newtonian fluids and their behaviors are out of interest. The solver efficiency was improved by decreasing the condition number of linear system using the change of variable.

Cracking is a challenge for processing ceramics. With the track irregularity, the simulation (Figure 4.13) is able to predict stress concentrations in the zones with smaller cross section. The advantage of an auxiliary laser in decreasing the maximum tensile stress and thus cracking risks is justified and evidenced by the numerical simulation (cases #2-4 in Section 4.5.2). Results also indicate that it is better to impose the auxiliary laser at the tail of melt pool where solidification begins. In addition, a laser with deeper penetration seems better than the one with limited irradiation on surface.

• Experimental comparison and other applications

In the comparison of melt pool shape, simulations show a coherence with experiments under different process conditions (Section 5.3.1). However, the comparison of track irregularity is very difficult. Despite this, the simulation can still reveal the tendency for certain cases (Section 5.3.2). In the application to multi-track deposition, it is shown that the hatch distance has significant influence on the layer surface (Section 5.4.2). Although the scanning directions may have also an impact, it is limited to the turning back region. By contrast, the zig-zag scanning leads to high temperature in the turning back region, meaning risks of potential evaporation. The preheating and reheating in multi-track deposition decreases the temperature gradient and consequently the maximum tensile stress.

Perspectives

This model provides the basic functionalities for the simulation of SLM process at track scale. However, there are some challenges during the development of this model and some of them still remain to be confronted or optimized. On the other hand, the model can be furtherly enriched and exploited. Hence, future work may be carried out in several directions as follows:

• Numerical resolution to be optimized The first issue is the condensation from powder to dense matter related to the assumption of continuous powder. In this model, an artificial condensation interval is assumed considering that instantaneous melting of powder induces high velocity and numerical convergence problems. The choice of this interval is delicate. The second issue is mass conservation. This needs a further investigation of the volume contraction rate θ in mass conservation, the level set method and mesh adaptation. A quick and efficient but not so elegant treatment is to compensate the mass variation by an artificial displacement of level set at gas/liquid interface. The third issue is the high computation charge of simulation. An optimization possibility can be addressed to the asynchronization of resolution, meaning to use different times steps for the solid state mechanical resolution on one hand and other resolutions on the other hand. Another possibility would be to reduce the mesh elements. For example, a body-fitted mesh [START_REF] Shakoor | A new body-fitted immersed volume method for the modeling of ductile fracture at the microscale: Analysis of void clusters and stress state effects on coalescence[END_REF] may be used for the non-impacted gas/material and powder/substrate interfaces.

• Application to other materials Application of this model to other materials is undergoing, in particular to metals, as they are largely used in industry. For that, material properties defining the heat source model should be adapted. One significant difference is that the viscosity of metals is usually one or two orders of magnitude lower than that of ceramics, meaning higher melt pool dynamics. Besides, the higher surface tension and Marangoni effect can also make contributions to higher melt pool dynamics. Consequently, attentions should be paid to the numerical resolution possibly requiring the use of turbulent models. In addition, processing metals is more prone to evaporation and keyhole due to their relatively lower boiling point. Hence, the integration of these phenomena is required.

• Other physical phenomena The evaporation and resulted recoil pressure and keyhole effect may be the first necessity to be integrated in the current model. Under this phenomenon, heat and mass are taken away by vapor forming at the liquid surface once attained the boiling point. The heat loss should be taken into account in the energy conservation and compressible behavior should be considered for the liquid-vapor phase change. This can be addressed by numerical treatment such as those proposed by Courtois et al. [START_REF] Courtois | A new approach to compute multi-reflections of laser beam in a keyhole for heat transfer and fluid flow modelling in laser welding[END_REF]. The recoil pressure is implicitly induced due to the expansion from liquid to vapor. Other contributions to melt pool dynamics like buoyancy and Darcy damping forces may be also taken into consideration.

An important issue in SLM process is the microstructure. In order to model this, the coupling of heat transfer with thermodynamic database (e.g. Thermocalc) is necessary. The modeling of microstructure by cellular automaton and finite elements (CAFE model [128]) was previously developed by Chen et al. [START_REF] Chen | Three-dimensional cellular automatonfinite element modeling of solidification grain structures for arc-welding processes[END_REF] in Cemef for welding. However, different to welding, the challenge of rapid solidification can be anticipated in the modeling of SLM. The integration of CAFE model will be a benefit to investigate the influence of process conditions (laser power, scanning speed and scanning strategy) on microstructral evolution.

• Coupling of two models at track and part scales Knowing that the developed model at track scale is difficult to be extended to the whole part, a suitable method to have this access is to couple it with macro models at part scale such as the one developed by Zhang et al. [START_REF] Zhang | Macroscopic thermal finite element modeling of additive metal manufacturing by selective laser melting process[END_REF] in Cemef. Before the coupling, the simplified heat source model in the macro model may be firstly validated by comparing with this model. For the coupling, one possible approach is to retrieve the stationary temperature field around the track from the micro model to the macro model. This may be carried out by extracting a specific subdomain (around the developing track) from the principal simulation with macro model. This domain is then used in the simulation with micro model to get precise temperature. The obtained temperature filed is representative and may be directly imposed in the macro simulation. Several difficulties may exist. The first is to determine when to launch the micro simulation, or equivalently when the macro simulation attains a stationary regime and when it is switched to the transient regime. The second difficulty is to choose the subdomain and the boundary conditions for the micro simulation. With this strategy, it is hopeful to have a full simulation at part scale with precise temperature evolution under a reasonable computation time.

R ésum é

L'application du proc éd é SLM est limit ée par la difficult é à contr ôler le proc éd é. Son application aux c éramiques est rendu particuli èrement difficile, en raison de la faible absorption de ces mat ériaux au rayonnement laser et de leur faible r ésistance au choc thermique. La maîtrise de ce proc éd é n écessite une compr éhension compl ète du transfert de chaleur, de la dynamique des fluides et de la m écanique des solides.

Abstract

The application of SLM process is limited by the difficulty of process control. Its application to ceramics is especially challengeable due to their weak absorption to laser and weak resistance to thermal shock. The mastery of this process requires a full understanding of heat transfer, fluid dynamics in melt pool and solid mechanics. In this work, we propose a numerical model for the simulation of SLM process applied to ceramics. The model is developed at the track scale and with the assumption of continuous powder bed. It is based on level set method and multiphase homogenization, with which we are able to follow the evolution of gas/material interface and phase transformation. Simulations are performed to study the influence of material properties and process parameters on temperature, melt pool shape, fluid dynamics and solid mechanics. Apart from the laser power and scanning speed, material absorption is also found to be important to the thermal behavior and the melt pool shape. With the fluid dynamics, convex shape of track cross section is achieved under surface tension. Besides that, liquid droplets collapsing formed by the melting of powder create melt pool instability when falling, thus leading to track irregularity after solidification. The Marangoni effect, caused by surface tension gradient at gas/material interface, is investigated. Its influence on temperature distribution, melt pool shape and track regularity is recognized. One interesting finding is the smoothing effect of track surface with negative ∂γ/∂T . When combine surface tension with scanning speed, track surface becomes more irregular with the increase of scanning speed. The well-known balling effect is reproduced with high scanning speed. This can be helpful to find the regime for regular track shape with given laser power and scanning speed. Cracking defect is deleterious in additive manufacturing. The use of an auxiliary laser can help to avoid this defect by decreasing the maximum tensile stress. The process mode of this auxiliary laser remains an interesting subject to be studied and some guidelines have been given by the presented simulations. The model is validated by the comparison of melt pool shape with experiments under different process conditions. Simulations can also reveal the tendency of track surface variation for certain cases. By the application to multi-track deposition, the influence of hatch distance on layer surface, temperature and stress evolution is emphasized.

1 . 1 1 . 2 1 . 3 1 . 4

 11121314 Additive Manufacturing . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Principle, advantages and challenges . . . . . . . . . . . . . . . . . . 1.1.2 Different processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . Selective Laser Melting . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Process principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Laser, materials and powder . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Physical phenomena in SLM . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Heat and mass transfer . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Melt pool dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Thermally induced stress . . . . . . . . . . . . . . . . . . . . . . . . 1.3.4 Microstructure and texture . . . . . . . . . . . . . . . . . . . . . . . Modeling of SLM process . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Micro modeling at track scale . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Macro modeling at part scale . . . . . . . . . . . . . . . . . . . . . . 1.4.3 Challenges . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Parts fabricated by different AM processes: (a) CLIP (Carbon3D [11]); (b) FDM [12]; (c) DMD (DMG MORI [13]); (d) SLM (similar for SLS) [14].

  .

  3 from step 1 to 3:

Figure 1 . 2 :

 12 Figure 1.2: (a) Inside view of SLM machine (SLM solutions [18]) and (b) 3D schematic drawing of the machine with name of main items [19].

Figure 1 . 3 :

 13 Figure 1.3: 2D schematics of the main processing steps for SLM. Step 1: After a layer of powder deposition; Step 2: construction; Step 3: ready for powder deposition. Green and blue color are used for powder with different apparent densities due to compaction during deposition and red color for condensed part.

Fe, Fe 3

 3 Al, Fe-Al intermetallics Fe-Ni, Fe-Ni-Cr, Fe-Ni-Cu-P 304L and 316L stainless steel H13 and H20 tool steel, M2 high speed steel Titanium based commercially pure Ti Ti-6Al-4V, Ti-6Al-7Nb, Ti-24Nb-4Zr-8Sn Nickel based Inconel 625 and 718, chromel, Hastelloy X Nimonic 263, IN738LC, MAR-M 247

1. 2 .Figure 1 . 4 :

 214 Figure 1.4: (a) New and (b) recycled nickel alloy powder in SLM [21].

Figure 1 . 5 :

 15 Figure 1.5: Parts fabricated by SLM: (a) jet engine [25]; (b) titanium orthopedic implant [26]; (c) 316L lattice structure for mechanical tests [27]; (d) 316L resistojet with channels [28].

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

 16117 Figure 1.6: Powder denudation around a single track in SLM with Ti-6Al-4V at ambient pressure varying from 0.5 to 220 Torr (left to right, 1 Torr = 133.3 Pa) [35].

Figure 1 . 8 :

 18 Figure 1.8: Cracks in SLM applied to Al 2 O 3 -ZrO 2 ceramics (provided by CdM, private communication): (a) top view, and (b) transversal cut view.

Figure 1 . 10 :

 110 Figure 1.10: Microstructures of Al 2 O 3 -ZrO 2 eutectic system by SLM provided by Centre des Matériaux. (a) fine eutectic lamellar, top view at track center, and (b) coarse eutectic lamellar, top view at the boundary between powder and solidified track. The process parameters are P L = 42 W and v L = 50 mm • s -1 .

Figure 1 .

 1 Figure 1.11: Grain orientation map obtained by electron backscattered diffraction (left) in front and (right) top view of an AlSi10Mg SLM part [47]. Dashed line indicates the boundaries of melt pool.

1. 4 .Figure 1 .

 41 Figure 1.12: Random powder bed: (a) schematic of the rain model for random packing with rotations; (b) powder bed produced by the rain model with higher relative density; (c) remove of some particles and non-physical holes; (d) cross section of a real powder bed (titanium alloy) [48].

Figure 2 . 1 :

 21 Figure 2.1: (a) Schema of two flux model in ray tracing; (b) Schema of Beer-Lambert law with only one forward flux. Red arrows for forward flux, blue arrows for backward flux and black arrows for the absorbed flux.

Figure 2 . 2 :

 22 Figure 2.2: Evolution of melt pool shape and temperature iso-contours. Red contour corresponds to the melt pool [49].

Figure 2 . 3 :

 23 Figure 2.3: Overhang case obtained by Hodge et al. [55]. The colors indicate the phase state. Blue corresponds to powder and red corresponds to consolidated material. The melt pool is shown with the black contour.

Figure 2 . 4 :

 24 Figure 2.4: Temperature histories for a point in the bulk (blue) and overhang region (red) [55].

Figure 2 . 5 :

 25 Figure 2.5: Comparison of melt pool cross section between experiment and different heat transfer models with scanning speed (a) 10 mm • s -1 and (b) 5 mm • s -1 and laser power of 1000 W. PLF -with latent heat of fusion and fluid flow, PHL -with latent heat of fusion, PH -pure heat conduction [56].

Figure 2 . 6 :

 26 Figure 2.6: Macro modeling at part scale developed by Zhang et al. [59] at Cemef: (a) decomposition of a layer by fractions for equivalent heating and cooling cycles; (b) simulation of heat transfer for the construction of a IN718 turbine with size of φ20 mm ×10 mm. The computation time is 17h with 60 Intel cores.

Figure 2 . 7 :

 27 Figure 2.7: Modeling schema of AM by SLM with the level set method to track gas/material boundary ψ = 0 [71].

Figure 2 . 9 :

 29 Figure 2.9: Optical indices of single Al 2 O 3 crystal, after [80].

Figure 2 . 10 : 3 L

 2103 Figure 2.10: Iso-contours of heat source distribution qL with different r int and α under the same R and P L . Iso-values correspond to 0.1, 0.2, 0.3 and 0.4 qLmax , where qLmax = (1 -R) 2P L πr 3 L

Figure 2 .

 2 Figure 2.11: One-way transformation from powder (zone Z 1 ) and dense matter (Z 2 ) with (a) volume fractions as a function of temperature in domain D 1 and (b) Illustration of possible distribution of zones and phases (same as structures) depending on the maximum heating temperature reached during thermal history.

Figure 2 .

 2 Figure 2.12: (a) Parallelization by FE mesh; (b) Parallelization by grid for 2D case. Red, blue and green represent the partition of FE mesh or grid points occupied by different cores. Black FE mesh or grid points are accessible for each core.

Figure 2 .

 2 Figure 2.13: 2D case: (a) Interpolation function N A (D) = S BCD S ABC from FE mesh to grid; (b) Integration in z direction based on grid; (c) Interpolation function N g G (M ) = S EP M Q

Figure 2 .

 2 Figure 2.14: Algorithm schema of computing qL , i, j ∈ [1, NbGrid], k, l ∈ [1, NbE].

Figure 2 .

 2 15(b).

Figure 2 .

 2 Figure 2.15: (a) Thermal conductivity of alumina λ a ; (b) Specific enthalpy of alumina h a .

2. 6 .

 6 Figure 2.16: Configuration of simulated system. The material domain contains only a dense alumina substrate and the rest is gas.

Figure 2 .

 2 Figure 2.17: (a) Heat source distribution qL in numerical simulation; (b) Comparison of qL along the laser axis between analytical and numerical solution.

Figure 2 .

 2 Figure 2.18: Iso-contours of qL from 0.5 to 4 (×10 14 W • m -3 ) with step of 0.5 × 10 14 W • m -3 .

4 :

 4 Test cases #1-4 to study the influence of material absorption. Red color indicates parameters changed with respect to reference case #1.

Figure 2 .

 2 Figure 2.19: System enthalpy evolution and energy received from laser in the case with adiabatic boundary conditions.

Figure 2 . 20 :Figure 2 . 21 :

 220221 Figure 2.20: Top view of temperature field evolution at the material surface in reference case #1 at t = 1, 2, 4, 6 ms. Iso-contours correspond to temperature T = 1000, 1500, 2004 (solidus), 2104 (liquidus, red) and 2500 • C.

Figure 2 . 22 :Figure 2 . 23 :

 222223 Figure 2.22: Longitudinal cut view of temperature distribution and temperature iso-contours corresponding to T = 1000, 1500, 2004 (solidus), 2104 (liquidus, red), 2500, 4000 and 6000 • C at t = 13 ms, for cases #1-4.

Figure 2 .Figure 2 .

 22 Figure 2.25: Longitudinal cut view of temperature distribution and temperature iso-contours corresponding to T = 1000, 1500, 2004 (solidus), 2104 (liquidus, red), 2500, 4000 and 6000 • C at t = 13 ms, for cases #5-7.

Contents 3 . 1

 31 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Modeling of fluid dynamics in melt pool . . . . . . . . . . . . . . 3.2.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Integrated driving forces . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Numerical resolution . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Semi-implicit formulation of surface tension . . . . . . . . . . . . . . 3.3.2 FE implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Remeshing strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Metric construction . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Intersection of metrics . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3 Blockage of elements . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Verification of energy and mass conservation . . . . . . . . . . . . 3.6 Single track simulation with powder on substrate . . . . . . . . . 3.6.1 Material properties and configuration . . . . . . . . . . . . . . . . . 3.6.2 Modification of convection velocity in heat transfer equation . . . . . 3.6.3 Results of reference case . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.4 Influence of material and process parameters . . . . . . . . . . . . . 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 :

 31 Figure 3.1: Track formation (a) with and (b) without surface tension applied to stainless steel with laser power 150 W and scanning velocity 5 m • s -1 . The powder layer thickness is 35 µm and it is deposited on a substrate with dimensions of 1 × 0.3 × 0.1 mm 3 . The color map varies linearly from room temperature (blue) to melt temperature (red) [49].

  (a) Layer 1, P L = 200 W, v L = 2 m • s -1 . (b) Layer 1, P L = 166 W, v L = 1.6 m • s -1 . (c) Layer 2, P L = 166 W, v L = 1.6 m • s -1 .

Figure 3 . 2 :

 32 Figure 3.2: Track stability at (a) high scanning speed; (b) low scanning speed (same linear energy) for first overhang layer and (c) low scanning speed on the (b) previous layer [69].

Figure 3 . 3 :

 33 Figure 3.3: Comparison of simulated (left) and experimental (right) melt pool shapes in a transverse cross section with P L = 150 W and v L = 400 mm • s -1 for TiC/AlSi10Mg material [44].

Figure 3 . 4 :

 34 Figure 3.4: 3D simulation of melt pool dynamics under (a) only constant surface tension; (b) Marangoni effect (surface tension included) and (c) Marangoni effect with recoil pressure. The color map corresponds to temperature between 293 (blue) and 4000 K (red). The red contour and arrows indicate the melt pool boundary and velocity direction, respectively [68].

Figure 3 . 5 :

 35 Figure 3.5: 3D simulation of pore formation illustrated in cross section: (a) depression of melt pool; (b) falling of liquid droplets and (c) close of connected melt and entrapment of gas [68].

Figure 3 . 6 :

 36 Figure 3.6: Convection flow induced by different Marangoni coefficients having opposite signs with T A and T C lower than T B (after [94]).

Figure 3 . 7 :

 37 Figure 3.7: Schematics of surface tension

Figure 3 . 8 :

 38 Figure 3.8: 2D geometrical interpretation of the metric tensor M i , which corresponds to an ellipse with semi-principal axis h m 1 and h m 2 along in direction m 1 and m 2 , respectively.

Figure 3 . 10 :

 310 Figure 3.10: Mesh adaptation in a 1 × 1 × 1 domain with respect to an ellipse (semi-principal axis 0.2 and 0.3) inclusion in a unit square. A Heaviside function is defined with a half thickness of 0.005. This function is used for the mesh adaptation with err = (a) 0.01; (b) 0.1; (c) 1; (d) 10; (e) 100; (f) 400. The imposed total element number is 10000 (the obtained number of elements is around 10000) and the minimum mesh size is set to 5 × 10 -4 .

  [START_REF] Chan | A two-dimensional transient model for convection in laser melted pool[END_REF] whereλ i = e T i M 1 e i and µ i = e T i M 2 e i (i= 1, 2, 3) and P = [e 1 e 2 e 3 ]. Then the intersection of 74 3.4. REMESHING STRATEGY two metrics are:

Figure 3 .

 3 Figure 3.12: Reinitialization of level set ψ: (a) ψ at t -∆t; (b) interface ψ = 0 deduced at t -∆t and distance of points A and E to ψ = 0; (c) ψ at t; (d) interface ψ = 0 deduced at t.

3. 5 .Figure 3 .

 53 Figure 3.13: Test case for energy and mass conservation: (a) system dimension for simulation with initial interface position (black line); (b) temperature and (c) velocity profile of test #2, at t=0, 2, 5, 10, 20, 30 ms.

Figure 3 .

 3 Figure 3.14: Conservation verification by the evolution of: (a) Energy (#2) and (b) mass.

  (a)) by their volume fractions.

Figure 3 . 15 :

 315 Figure 3.15: Physical properties of material with (a) conductivity of alumina powder λ Z 1 with 50% porosity and (b) viscosity of powder η Z 1 and dense matter η Z 2 .

Figure 3 . 16 :

 316 Figure 3.16: Temperature field snapshot (left) and profile (right) along the white line at the beginning of laser scanning: (a) without and (b) with the modification of convection velocity using Eq.3.54 in heat transfer.

Figure 3 . 17 :

 317 Figure 3.17: Snapshots of reference case #1 at t = 5 ms with laser at position X = 1.2 mm: (a) temperature and iso-contours for T = 300, 700, 1100, 1500, 1900, 2104, 2300, 2700 and 3100 • C; (b) mesh around the solidified track at position near X = 0.2 mm; (c) mesh around melt pool near laser position; (d) heat source distribution qL and different contributions to mesh adaptation. Top view without gas for (a-c) and longitudinal cut view for (d).

Figure 3 .

 3 Figure 3.18: Blockage of elements in different regions: region A -solidified gas/material interface, region B -non impacted (by laser) gas/material or powder/substrate interface. The magenta element is not changed at t = 6, 8 and 10 ms.

Figure 3 .

 3 Figure 3.19: Comparison of solidified track surface: (a) with and (b) without blockage at t = 2, 3, 3.6 ms. The red iso-contour corresponds to T l = 2104 • C.

Figure 3 . 20 :

 320 Figure 3.20: Formation of droplets during t = 6.08 to 6.24 ms. Arrows indicate the velocity direction and magnitude by color truncated to 500 mm • s -1 . The red contour corresponds to T l = 2104 • C and the black one to the region exposed to laser irradiation.

Figure 3 . 21 :

 321 Figure 3.21: Verification of conservation of (a) mass and (b) energy for reference case #1.

3. 6 .

 6 SINGLE TRACK SIMULATION WITH POWDER ON SUBSTRATE Top view of case #2-5.

  Transversal cut view of case #2 and #3.

Figure 3 . 22 : 89 CHAPTER 3 . 3 . 6 . 4 . 2 Figure 3 . 23 :

 3228933642323 Figure 3.22: Melt pool (red contour) and liquid velocity field when the laser (black contour) arrives to X = 1.416 mm in (a) top and (b) transversal cut view. White line indicates the initial powder/substrate interface.

Figure 3 . 24 :

 324 Figure 3.24: Track height H app of case #1-5 with respect to the substrate surface. Isocontours correspond to ∆z = 15, 20, 22.5, 25, 27.5, 30 and 35 µm. The height variation can be estimated by the separated iso-contours. Note that it is scaled by a factor 2 only in Y direction for better visualization.

Figure 4 . 1 :

 41 Figure 4.1: Stress distribution: (a) σ xx and (b) σ zz at the end of the 16 th layer [116].

Figure 4

 4 .2; • temperature in the remelted zone is considered to be uniform and equal to melting point and it equals to the ambient temperature in other zone.

Figure 4 . 2 :

 42 Figure 4.2: Stress field around the remelted zone in stable regime with different initial temperatures when processing Al 2 O 3 [117]. (top) T 0 = 0 • C; (bottom) T 0 = 1600 • C. Note that the shape of remelted zone is presupposed.

Figure 4 . 3 :

 43 Figure 4.3: Stress distribution in scanning direction compared with nominal process conditions of P L = 840 W, r L = 200 µm, v L = 0.5 m • s -1 and T 0 = 650 • C when processing Ti -6Al -4V: (a) different laser sizes scaled by factor ω f and power scaled by P f while keeping the same maximum surface laser intensity; (b) laser power increased by 20%; (c) scanning speed decreased to 0.1 m • s -1 ; (d) different initial temperatures [111].

Figure 4 . 4 :

 44 Figure 4.4: (a) 2D element E decomposed to sub-elements E i , E j and E k ; (b) Bubble function N b (x) in the element E; (c) Bubble function in sub-element E i related to linear interpolation function by N b (x) = DN i (x).

Figure 4 .

 4 4(a). The bubble function has a value of 1 at G and vanishes at element boundary ∂Ω E as shown in Figure 4.4(b). By isolating the sub-element E i , one can relate the bubble function N b (x) in this sub-element to linear interpolation function by N b (x) = DN i (x) as shown in Figure 4.4(c). Another important property of bubble function is:

  27) where the notations ( * : * ) = Ω h * : * dΩ, ( * , * ) = Ω h * • * dΩ and ( * , * ) Γ h f = Γ h f * • * dΓ are applied to simplify the expression. The first equation can be separated by isolating the linear 105 CHAPTER 4. SOLID MECHANICS test function and the bubble test function:

4 . 39 B 0

 4390 Eq.4.42, σ Eq.4.11 B 0 < σ s = s elP , ¯ = ¯ t- ˙ Eq.4.41 (Newton Raphson) σ Eq.4.11 ⇒ λ Eq.4.12 ⇒ s Eq.4.39 C, C b Eq.4.43,4.45 T ⇒ K, m s Eq.4.15

Figure 4 .

 4 Figure 4.6 shows global flow chart of resolution procedure.Firstly, it should be pointed out that at each time step, all the heat transfer, fluid dynamics, transport and reinitialization of LS, solid mechanics and update of coordinates are solved successively. Remeshing is executed with a fixed frequency but not at each time step. The frequency is chosen in the way that the gas/material interface does not exceed the refined transition zone in order to assure good interface representation. Thus, it depends on the time step, the velocity of interface and the size of refined transition zone. Note that it may be adapted according to different melt pool dynamics, for example, when processing metals with much lower viscosity.The heat source qL should be firstly obtained as an input for heat transfer. This part is detailed in Section 2.5.1. Together with the convection velocity u - f and temperature T - of previous time step as input, the current temperature distribution can be obtained via the temperature solver (Section 2.5.2). Phase fractions and material properties are then calculated according to temperature. They are used to evaluate the volume variation rate θ

Figure 4 . 6 :

 46 Figure 4.6: Algorithm flow chart of global resolution (after [3]) coupling heat transfer, fluid dynamics and solid mechanics (Figure 4.5).

Figure 4 Figure 4 . 7 :

 447 Figure 4.7: Simulated domain with mesh for static welding.

Figure 4 .Figure 4 . 8 :

 448 Figure 4.8: Equivalent stress at 10 s: (a) in symmetric plane (gas removed); (b) along radial profile (red line in (a)) of σ at metal surface.

Figure 4 . 9 :

 49 Figure 4.9: EVP behavior assumed for Al 2 O 3 ceramic: (a) stress-strain relation at T = 20, 200, 400, 600, 800, 1000 and 1200 • C by interpolating the results of Sánchez-González et al. [124] and assuming m = 0; (b) sensibility coefficient assumed to have same tendency as in [3].

Figure 4 .Figure 4 .

 44 Figure 4.11: The mesh blocked in remelted zone behind melt pool at t = 2, 3 and 4 ms in case #1. Red contours correspond to T = 2004 and 1804 • C (displayed only in material domain).

Figure 4 .

 4 Figure 4.12: Stress distribution σ xx , σ yy and σ zz at t = 8 ms (laser at X = 2.6 mm) in case #1. The black and horizontal white lines indicate the gas/material and initial powder/substrate interfaces, respectively. White iso-contours correspond to T = 500, 1000, 1500 and 2004 • C.

Figure 4 .

 4 Figure 4.13: Stress distribution σ xx , σ yy and σ zz , track height (right y axis) along profile (horizontal white line in Figure4.12) and track surface (with powder) in case #1. Note that the temperature increases along X and the laser is at X = 2.6 mm at this time t = 8 ms.

Figure 4 .

 4 Figure 4.14: Stress distribution σ xx , σ yy and σ zz along a vertical profile in the longitudinal symmetric plane (vertical red line in Figure 4.12) in case #1. The bottom of remelted region is indicated by the vertical dashed magenta line. Vertical black line is the position of initial powder/substrate interface.

Figure 4 . 15 :Figure 4 . 16 :

 415416 Figure 4.15: Melt pool shape of cases #1-4 (black, red, blue and green, respectively) in longitudinal symmetric plane at t = 8 ms.

Figure 4 . 17 :

 417 Figure 4.17: Stress distribution σ xx in case #3 (2nd Yb:YAG laser behind 1st laser by 250 µm) at time t = 8 and 12 ms. Note that there is no heat input after t = 8.667 ms.

Figure 4 .

 4 Figure 4.18: Stress distribution σ xx in case #4 (2nd CO 2 laser behind 1st laser by 250 µm) at time t = 8 and 12 ms. Note that no heat input after t = 8.667 ms.

Figure 4 .

 4 Figure 4.19: Stress distribution σ xx in cases #1-4 along red profiles indicated by vertical red lines in Figure 4.12, 4.16 and 4.17-4.18. The vertical dashed line shows the bottom of remelted zone (red color for case #2 and blue color for cases #1, 3 and 4). Vertical black line is the position of initial powder/substrate interface.

  Figure 4.19: Stress distribution σ xx in cases #1-4 along red profiles indicated by vertical red lines in Figure 4.12, 4.16 and 4.17-4.18. The vertical dashed line shows the bottom of remelted zone (red color for case #2 and blue color for cases #1, 3 and 4). Vertical black line is the position of initial powder/substrate interface.

Figure 5 . 1 :

 51 Figure 5.1: Crushed Al 2 O 3 powder particles with 1wt% of carbon absorbers heterogeneously distributed on particle surfaces (CdM). Powder particles are fully or partially covered by absorbers and some of them are not covered at all.

Figure 5 . 2 :

 52 Figure 5.2: Schematic diagram of single track scanning of a layer of Al 2 O 3 powder deposited on Al 2 O 3 substrate for laser power P L = 42, 84, 126 and 168 W. For each power, the scanning speed is varied in 50, 100, 200, 400 and 800 mm • s -1 .

Figure 5 . 3 : 129 CHAPTER 5 .Figure 5 . 4 :

 53129554 Figure 5.3: Height colormap of track and powder after single pass scanning with different velocities and same laser power P L = 84 W. The height increases from green to yellow color.

Figure 5 . 5 :

 55 Figure 5.5: Transverse sections observed in SEM.The transverse melt pool shape corresponds to the boundary of melted zone which can be noticed as it is stripped from the non-melted zone. Besides the 5 cases at the bottom left with different highlighted scales, a scale of 100 µm is used in all other cases. Melt pool shape of cases in red rectangles with linear energy E l between 0.2 and 0.9 J • mm -1 is considered more suitable for the construction. Note that the melt pool shape for a given condition may vary with different cutting positions.

Figure 5 . 6 :

 56 Figure 5.6: Definition of melt pool width w app , depth H zr and height H app .

Figure 5 . 7 :

 57 Figure 5.7: Experimental results of melt pool dimension evolution with respect to scanning speed with different laser power: (a) width and (b) depth. The measurement uncertainty is 5 µm, systematically.

Figure 5 . 8 :

 58 Figure 5.8: Integrating sphere for measuring reflection: (a) instrument and (b) schematic diagram of uniform scattering on inner surface.

Figure 5 . 9 :

 59 Figure 5.9: Experimental result of reflection measurement and interpolation by a power law.

5. 2 .Figure 5 .

 25 Figure 5.11: Comparison of melt pool (a) width and (b) depth between experiments (square marker with error bar), analytical prediction (dashed line with cross markers) and simulation (solid line with circle marker) with calibrated heat source parameters. Black, red and blue colors are used for P L = 42, 84 and 126 W, respectively.

  )

  Figure 5.12: Comparison of melt pool shape for different process conditions. The experimental results are overlaid with simulated melt pool contours (red profiles) with the same scale. Note that the melt pool shape for a given condition may vary with different cutting positions.

Figure 5 .

 5 Figure 5.13: Comparison of melt pool (a) width and (b) depth between experiments (square marker with error bar) and simulation (solid line with circle marker) with optimized heat source parameters calibrated by Eq.5.18. Black, red and blue colors are used for P L = 42, 84 and 126 W, respectively.

Figure 5 .

 5 Figure 5.14: Comparison of track height map and profile (in longitudinal median plane) between experiments and simulations for cases (a) P42V200 and (b) P84V400. Only the stable regime of height profile is plotted for simulation results.

Figure 5 . 15 :

 515 Figure 5.15: Balling effect observed in simulation with v L = 600 mm • s -1 and P L = 84 W.

Figure 5 . 16 :

 516 Figure 5.16: Balling in processing 316L with SLM ([36]) and simulated temperature distribution at different instants for the case P84V600. The black contour corresponds to the melt pool. The graph at the bottom shows the temperature profile at t = 4.4 ms and Z = 0.94 mm (9 µm below the initial substrate surface) in the median XZ plane.

Figure 5 . 17 :

 517 Figure5.17: Configuration of simulated system. The material domain contains a dense alumina substrate and a layer of powder. The rest is gas.

Table 5 . 2 :

 52 Figure 4.6, including the update of mesh nodes with the velocity obtained by the resolution of solid mechanics. The problem of element distortion mentioned in Section 4.5 is avoided by furtherly increasing the CHAPTER 5. EXPERIMENTAL VALIDATION AND OTHER APPLICATIONS Case Different scanning strategies in multi-track deposition for cases #I-III.

9 YFigure 5 .

 95 Figure 5.18: Schematic diagram of multi-track deposition with unidirectional scanning andP L = 84 W, v L = 400 mm • s -1 .The hatch distance ∆y is 50 µm for case #I, II and 40 µm for case #III. Variable evolution at 9 points (black, red and blue) will be investigated. They are located at Z = 0.94 mm (9 µm below the initial substrate surface) and numbered in circle from 1 to 9 from pass #1 to #3 and from left to right. Pass #1 is always at y = 0.25 mm and pass #2 is in inverse direction for zig-zag strategy (case #II).

#Figure 5 .

 5 Figure 5.19: Colormap of layer height at the end of (a) pass #1 (same for all cases); (b) pass #2 and (c) pass #3 for cases #I-III with laser condition P84V400.

Figure 5 . 20 :

 520 Figure 5.20: Profile of gas/material interface in cross section (X = 1.5 mm) at the end of each pass for cases (a) #I; (b) #II and (c) #III. Vertical dashed lines indicate the Y position of laser in each pass.

  #III, unidirectional, ∆y = 40 µm

Figure 5 . 22 :

 522 Figure 5.22: Temperature distribution at different instants for cases (a) #I-III; (b) #I; (c) #II and (d) #III. The white line for cases #I and III at t = 9.75 corresponds to the scanning line (Y = 0.25 mm) of pass #1. White arrows with anchor point indicate the laser position and moving direction at different instants.

Figure 5 . 23 :

 523 Figure 5.23: Axial stress distribution σ xx in the XZ plane (Y = 0.25 mm, the scanning line of pass #1) at different time for cases (a) #I-III; (b) #I and (c) #II. The red contour indicates the melt pool boundary and the white line corresponds to the initial substrate surface. Black arrows with anchor point indicate the laser position and moving direction at different instants.

Figure 5 . 24 :

 524 Figure 5.24: Temperature history at points 1-9 in case #I. The numbering of points and theirs locations are illustrated in Figure 5.18 with corresponding colors.

Figure 5 .

 5 Figure 5.25: Temperature history at points 1-9 in case #II. The numbering of points and theirs locations are illustrated in Figure 5.18 with corresponding colors.

Figure 5 Figure 5 .

 55 Figure 5.26: σ xx history at points 1-9 in case #I. The numbering of points and theirs locations are illustrated in Figure 5.18 with corresponding colors.

Figure 5 . 27 :

 527 Figure 5.27: σ yy history at points 1-3 in cases: (a) #I and (b) #II. The numbering of points and theirs locations are illustrated in Figure 5.18 with corresponding colors.

Figure 5 .Figure 5 . 28 :

 5528 Figure 5.28: σ xx history at points 1-9 in case #II. The numbering of points and theirs locations are illustrated in Figure 5.18 with corresponding colors.

  D i , • Z i , • S iHomogenized property in domain, zone and structure

	List of Symbols		List of Symbols
	Other Symbols I	Identity tensor
	K • f + s		Stiffness matrix =f if f > 0 and =0 otherwise Deviatoric stress tensor	Pa
	{•} e i , i = 1,2,3		Global property Basis
	( * , * ) f M	2,3	Eigenvalues = Ω ( * • * ) dΩ Marangoni force	Pa
	λ L ( * : * ) f s		Laser wavelength = Ω ( * : * ) dΩ Surface tension force	m Pa
	µ ∆ s f v		Shear modulus Surface Laplace operator Volumetric force	Pa N • m -3
	ν erf (x) f imp		Poisson's ratio Error function Imposed surface force	Pa
	Ω ω Ω E Ω φ B h H P h T h V h ∇ s tr( * ) G R		Simulated domain Angular frequency Domain occupied by an element Laser flux intensity Bubble velocity space Heaviside function Pressure space Mesh topology Velocity space Surface Nabla operator Trace of tensor Residual	s -1 W • m -2 m • s -1 Pa m • s -1
	ψ Roman Symbols u	Level set function Velocity field	m m • s -1
	ψ l n u imp		Level set function associated to melt pool boundary Optical index Imposed velocity	m m • s -1
	ρ ∆h v		Density Variation of volumetric enthalpy	kg • m -3 J • mm -3
	σ r ∆h v,f		Stefan-Boltzmann constant Volumetric enthalpy required for the fusion	W • m -2 • K -4 J • mm -3
	σ Y qL		Nominal yield stress Heat source distribution	Pa W • m -3
	τ K , τ C qr		Stabilization coefficients Heat lost rate by radiation	W • m -3
	C		Strain tensor 4th order tangent modulus
	σ E		Stress tensor Stiffness tensor	Pa Pa
	ε M		Emissivity Metric tensor	m -2
	P		Coefficient for change of variable Tensor composed by basis
	ϑ X		Refractive index Length distribution tensor	m 2
	Ω ė		Directional solid angle Deviatoric of strain rate tensor	s -1
			ix x

h Discretized simulated domain by FE ∂Ω f , ∂Ω u Boundary with imposed force or velocity ∂Ω T , ∂Ω q , ∂Ω cr Boundary with imposed temperature, flux or convection and radiation i Gradient at node i G ij = G j -G i Gradient difference between node j and i m i , i = 1,2,3 Principal orthogonal basis n Unit normal vector at gas/material interface n Ω Unit normal vector on system boundary n l Unit normal vector on melt pool boundary

  [START_REF] Li | Comparison of volumetric and surface heating sources in the modeling of laser melting of ceramic materials[END_REF], applied to Al 2 O 3 -based refractory during CO 2 laser treatment. They confirm experimentally that the volumetric heat source is more accurate than surface heat source in the prediction of melt pool shape for Al 2 O 3 -based ceramics. The same model is also used by Defillon et al.[START_REF] Defillon | Modèle du bain de fusion de SLM sur céramique avec disparition de l'agent absorbant restreinte au bain[END_REF] for Al 2 O 3 -ZrO 2 ceramics, with a Yb:YAG laser. Considering the assumption of continuous powder and the simplicity of this model, it will be adapted and used in our modeling hereafter. More details will be provided in section 2.3.2.

			2.1. STATE OF THE ART
	absorption coefficient α:
	dφ dz	= -αφ	(2.3)
	Korner et al. [48] used this model in the modeling at powder scale with Lattice Boltzmann
	method. In order to take into account the porous effect of powder, they introduce the material
	fraction g m to have the apparent absorption of powder αg m , leading to:
	∆φ ∆x	= -αg m φ	(2.4)
	where ∆φ is the energy absorbed within a cell size ∆x. A continuous model is used by Li
	et al.
				2)
	where P Ω → Ω is the probability function of a flux propagating in direction Ω scattered
	to direction Ω. One can note that this equation takes into account the scattering in all
	directions. Eq.2.2 can be simplified to Eq.2.1 if only one scattering direction is considered.
	Their results with metallic powder show the coherence with ray tracing simulations. A good
	agreement with experiments is achieved in the correlation between effective absorptance of
	powder and intrinsic absorptance of corresponding dense metals. Their model is used by
	Hodge et al. [55] and King et al. [69].
	Eq.2.1 can be furtherly simplified by ignoring the backward flux and scattering, leading
	to the Beer-Lambert law, which takes into account the local absorption of material. In this
	situation, the flux intensity decreases logarithmically when it penetrates into a medium with

Relocalization Search containing element (s) Other (s) Total (s)

	Yes	2.287	0.322	2.609
	No	∼ 0	0.148	0.148

Table 2 .

 2 3: Computation time for the cases with and without relocalization. Simulation with NbE = 397611 and NbGrid = 317520 at 28 Intel cores.

Table 3 .

 3

		# l (µm)	(µm)
		1	50	250
		2	30	250
		3	20	250
		4	20	150
		5	20	100
	1: Material properties	Table 3.2: Test Parameters

Table 3 .

 3

	1	-8.2	Eq.3.51	200
	2	-48	Eq.3.51	200
	3	48	Eq.3.51	200
	4	-8.2	0.2	200
	5	-8.2	Eq.3.51	400

3: Values of material and process parameters to study the melt pool dynamics with red color indicating changes with respect to reference case #1. Note that the viscosity obtained by Eq.3.51 at 2104 • C is about 0.04 Pa • s.

  .22.

		Unit	#1	#2	#3	#4	#5
	L	µm	35	40	31	36	20

Table 3 .

 3 4: Marangoni numbers for case #1-5 calculated by taking L as the half width of melt pool and ∆T = T max -T l . The parameters used to calculate the thermal diffusivity are λ

	8.2	-42	42	-8.2	-8.2

Table 4 .

 4 1: Number of iterations (resolution of linear system) of first Newton-Raphson iteration at different time steps and corresponding computation time with 24 Intel cores.

	Time step	Nb of iterations	Time ( s)
	(∆t)	Without	With	Without	With
	1	1334	1338	4.925	4.809
	2	1440	110	5.212	0.402
	50	1498	40	5.361	0.152
	100	1089	23	3.930	0.089
	150	1176	22	4.213	0.085
	200	968	20	3.476	0.078
	Total			1530	57

Table 4 .

 4 3: Total computation time for a simulation with 4333 ∆t until the extinction of principal laser and the time percentage consumed by each part. Note that the number of elements evolves from about 1 200 000 to 1 780 000. The calculation is carried out with 48 Intel cores.

	Total Heat source Heat transfer Fluid mechanics Solid mechanics Remesh Other
	(h)			(%)
	78.7	3.6	10.4	20.1	23.2	24.9	17.9

Table 5 .

 5 1: Process conditions used in set G for the calibration of heat source parameters.

	P L ( W)	42	42	42	84	84	84	126	126
	v L mm • s -1	50	100 200 100 200 400 200	400
	E l = P L /v L J • mm -1	0.84 0.42 0.21 0.84 0.42 0.21 0.63 0.315

  Dans ce travail, nous proposons un mod èle num érique pour la simulation du proc éd é SLM appliqu é aux c éramiques. Le mod èle est d évelopp é à l' échelle du cordon et avec l'hypoth èse d'un lit de poudre continu. Il est bas é sur la m éthode level set et l'homog én éisation multiphasique, avec laquelle nous sommes capables de suivre l' évolution de l'interface gaz/mati ère et les transformations de phase. La simulation d évelop ée permet d' étudier l'influence des propri ét és du mat ériau et des param ètres du proc éd é sur la temp érature, la forme du bain liquide, la dynamique des fluides et la m écanique des solides. En dehors de la puissance du laser et de la vitesse de balayage, l'absorption du mat ériau est également importante pour la thermique et la forme du bain liquide. En int égrant la dynamique des fluides, la forme convexe du cordon est obtenue en consid érant la tension de surface. Les gouttelettes liquides se forment lors de la fusion de la poudre et cr éent une instabilit é du bain. Ceci entraîne une irr égularit é du cordon apr ès solidification. L'effet Marangoni, provoqu é par le gradient surfacique de la tension de surface, est étudi é. Son influence sur la r épartition de la temp érature, la forme du bain liquide et la r égularit é du cordon est évoqu ée. Cet effet peut lisser la surface du cordon avec ∂γ/∂T n égatif. En augmentant la vitesse de balayage, la surface du cordon devient plus irr éguli ère. L'effet de "balling" est reproduit avec une vitesse de balayage élev ée. Cela peut être utile pour trouver le r égime donnant une forme de cordon r éguli ère étant donn ées la puissance et la vitesse du laser. Le d éfaut de fissuration est d él ét ère dans la fabrication additive. L'utilisation d'un laser auxiliaire peut aider à éviter ce d éfaut en diminuant la contrainte de traction maximale. Le mode de fonctionnement de ce laser auxiliaire reste un sujet int éressant à étudier et quelques pistes ont ét é donn ées par les simulations pr ésent ées. Le mod èle est valid é par la comparaison de la forme du bain liquide avec des exp ériences dans diff érentes conditions de proc éd é. Les simulations peuvent également r év éler la tendance de variation de la surface du cordon dans certains cas. En simulant le d ép ôt cons écutif de diff érents cordons, l'influence de taux de recouvrement sur la surface d'une couche, le champ de temp érature et l' évolution des contraintes est soulign ée.

	Mots Cl és
	Fabrication additive, mod élisation, level
	set, él éments finis, dynamique des fluides,
	m écanique des solides

Acknowledgement

4.10. The spot size of auxiliary laser is augmented to cover a larger region and the power is chosen such that the maximum temperature with only auxiliary laser attains around 1600 • C as in [START_REF] Hagedorn | Net shaped high performance oxide ceramic parts by selective laser melting[END_REF]. The auxiliary laser type will be also studied as Al 2 O 3 has higher absorption coefficient with CO 2 laser than Yb:YAG laser (Figure 2.9 [START_REF] Faure | Etude de l'interaction Rayonnement-Matière Dans Un Milieu Granulaire En Vue de l'application Au Procédé de Frittage Laser[END_REF]), leading to different heat inputs. Hence, in case #4 the absorption is increased while the power is decreased in order to keep almost the same maximum temperature with only auxiliary laser.

The system configuration is the same as in the previous chapter. The bottom is fixed and other surfaces are free for boundary conditions of solid mechanics [START_REF] Mukherjee | An improved prediction of residual stresses and distortion in additive manufacturing[END_REF]. The time step is still taken as ∆t = 2 µs. All simulations follow the resolution flow chart in Figure 4.6 except the deactivated coordinate update with velocity resulting from the resolution of solid mechanics.

Chapter 5

Experimental validation and other applications Contents

MEASUREMENT AND CALIBRATION OF PARAMETERS DEFINING THE HEAT SOURCE

(see the next section). This assumption have only slight influence as the effective input power (1 -R)P L does not change a lot. In order to take into account the influence of linear energy on the reflection, a power law is supposed for both powder and liquid, expressed as follows:

This model will be used in the following analytical model and numerical simulations.

Analytical model for prediction of melt pool shape

Here is presented an analytical model for the prediction of melt pool shape. This model was proposed by Defillon et al. [START_REF] Defillon | Modèle du bain de fusion de SLM sur céramique avec disparition de l'agent absorbant restreinte au bain[END_REF] and will be used to make a first calibration of parameters r int and α. In this model, some assumptions are made as follows:

• The laser scans unidirectionally a dense flat substrate in direction x with fixed power P L and speed v L . • There is no heat exchange by conduction or convection between a local point and its neighbors. • Absorption coefficient is distinguished between liquid and solid material, noted as α l and α s (Figure 5.10), respectively. Under these assumptions and using the Beer-Lambert absorption law, the volume heat source term ( W • m -3 ) in Eq.2.22 can be expressed by:

Scanning direction x

where Z l (x, y, t) describes the boundary of melt pool to be deduced as shown in Figure 5.10 and φ 0 is still the Gaussian distribution (Eq.2.21) but expressed in the laboratory referential

Appendix

The integration of the heat source distribution qL (second case inEq.5.3) from t 1 to t 2 is:

The integration term can be calculated by the change of variable:

where the error function erf (x) is defined by: erf (x) = 2 √ π

x 0 e -u 2 du (A.3) with erf (0) = 0, erf (-∞) = -1 and erf (+∞) = 1. Considering that the laser scans from t 1 = 0 and to t 2 → +∞, the energy absorbed by a point in the middle of the laser trajectory (x = v L t 2 /2) can be calculated by:

, y, 0, t 2 dt = lim