Finite element modeling of deposition of ceramic material during SLM additive manufacturing

. A three dimensional model for material deposition in Selective Laser Melting (SLM) with application to Al (cid:2) O (cid:3) - ZrO (cid:2) eutectic ceramic is presented. As the material is transparent to laser, dopants are added to increase the heat absorption efficiency. Based on Beer-Lambert law, a volumetric heat source model taking into account the material absorption is derived. The Level Set method with multiphase homogenization is used to track the shape of deposed bead and the thermodynamic is coupled to calculate the melting-solidification path. The shrinkage during consolidation from powder to compact medium is modeled by a compressible Newtonian constitutive law. A semi-implicit formulation of surface tension is used, which permits a stable resolution to capture the gas-liquid interface. The formation of droplets is obtained and slight waves of melt pool are observed. The influence of different process parameters on temperature distribution, melt pool profiles and bead shapes is discussed.


Introduction
The Selective Laser Melting (SLM) process is a technique of additive manufacturing which heats and melts powder layer by layer with an energy concentrated laser beam. During each procedure, the feed piston moves upward a certain distance and the build piston moves downward a layer thickness (Figure 1). Then the roller deposes a new layer of powder and laser beam begins to scan the powder bed according to a trajectory predefined by a 3D model. The laser melts the powder which solidifies to form a new solid material layer constituted of multiple beads. The procedure 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 density (up to 100%, Hagedorn et al. [1]) can be obtained with SLM, which means better mechanical properties. However, energy is so focused that process control becomes difficult. High thermal gradients usually exist and may result in defects such as cracking.
Ceramic material with eutectic composition of Al ଶ O ଷ − ZrO ଶ system (Figure 4) is studied. This material is potentially attractive for diverse applications (aeronautics, medical). It is light, refractory, and shows high resistance to creep. Regarding SLM processing, unlike metallic materials, which absorb most part of nonreflected laser energy on powder surface, it is transparent to Nd:YAG laser with wavelength 1064 nm. Dopants such as carbon particles are added in powder in order to increase and control absorption. Therefore, the heat source is on volume rather than on surface. Thermal phenomena in SLM are more complicated than in other processes. Unlike other traditional metallurgy processes, melting and consolidation take place locally and rapidly in SLM. A high cooling rate of 10 K/s was reported by Salehi and Dehghani [3] and refined phase size was obtained. Chou et al. [4] used pulsed laser beam rather than continuous one to increase the cooling rate. This resulted in higher undercooling and a significant microstructure deviation from equilibrium. The final part showed an increase in hardness compared with traditional cast alloy and other SLM studies.
Concerning the modeling of SLM processing, most researches are interested in heat transfer and prediction of bead shape. Hodge et al. [5] chose a complicated heat source model described by Gusarov et al. [6]. The temperature field was predicted, as well as the shape of melt pool. An interesting overhang case (melt pool over the powder rather than the constructed part or the substrate) was demonstrated, where two separate pools (in the overhung and non-overhung region) were obtained. They claimed that this phenomenon was caused by the insulating behavior of the unconsolidated powder under the overhung region. Li and Gu [7] predicted the temperature field in pure Ti − 6Al − 4V powder bed, with a Gaussian heat source model and temperature dependent thermal properties. They found that most heat was eliminated through conduction in the cold substrate and that the material experienced a rapid quenching process. Roberts et al. [8] used the heat source model considering the absorption of laser energy described in [9] and employed the method of element birth and death to present the addition of multiple layers, also for Ti − 6Al − 4V material. The same method is used by Marion et al. [10]. Quenching effect is also obtained in both [8] and [10].
In the present work, we aim to predict the temperature distribution, the bead formation and the dimension of melt pool. Their sensitivity to process parameter (scanning velocity) and material property (conductivity) is also studied.

Modeling method
In the present work the modeling of SLM processing is at the scale of the development of a single bead, for which a schematic drawing is illustrated in Figure 2. The scale is over the powder particle scale (several micrometers), so that continuous medium is considered for the material during the process. The modeling method should be able to track the bead shape during melting and solidification. The thermal resolution should take into account the phase transformation and use a dedicated heat source model, accounting for the interaction between laser radiation and the ceramic material in order to predict the temperature distribution.

Level Set method
Level Set (LS) method [11] is used to model the SLM processing. As shown in Figure 3, we consider the substrate (or, alternatively the previously consolidated solid material) and the powder bed as a material domain ‫ܦ(‬ ). It is surrounded by the gas domain ‫ܦ(‬ ଵ ). These two domains are separated by an interface ߰ = 0, where ߰ = ‫,ܠ(߰‬ ‫)ݐ‬ is the signed distance (to the interface) function of position ‫ܠ‬ and time ‫.ݐ‬ Centered on this interface, we define a transition zone with half thickness ߳, in which property (density, conductivity, etc.) of each domain is averaged by a Heaviside function ℋ to get global property ߯ [12]: where 〈߯〉 (݅ = 0,1) is the multiphase average property in domain ‫ܦ‬ and ℋ is defined by: The choice of ߳ is delicate because it should be small enough to capture precisely the interface and also big enough to have a smooth transition of property regarding the numerical resolution. The distance function ψ should be adapted at each time step in order to track the interface. This step is just after the resolution of Navier-Stokes (NS) equation. It corresponds to the solution of the following equation: where ‫ܝ‬ ௌ denotes the velocity of the interface. However, this process does not guarantee the equality between the geometrical distance ݀ from a point (except points at interface) to the interface and the value of |߰|. Therefore, ߰ is recalculated by a geometrical method [13] with respect to the new position of the interface resulting from Eq.(3), considering that the resolution of Eq.(3) is done only for the transportation of LS interface (߰ = 0).

Multiphase modeling
The material domain ‫ܦ(‬ ) is a multiphase system as several phases may coexist according to the phase diagram in Figure 4. Figure 5 shows the melting and solidification steps in material domain. At initial state ܶ < ܶ = 1860 °C, there is only powder with 3 phases in ‫ܦ‬ (step 1). Once temperature attains the melting temperature ܶ , powder begins to be melted and liquid fulfills gaps between particles. At this stage, ܶ keeps constant as the system is in eutectic composition and latent heat is released (step 2). Once powder is totally melted (step 3) and begins to solidify, solid material with eutectic composition is obtained (step 4). During SLM, those four steps are chained extremely rapidly, with possible partial melting and solidification of powder.  Following the same procedure as in [15][16], we define zones, structures and phases in material domain as shown in Table 1 and Table 2.
It should be pointed out that the transformation from Z to Z ଵ is non-reversible as the powder can never be recovered once melted. On the other hand, in Gas domain (D1), there is only one zone, one structure and one phase (g). They are not listed in Table 1.

Domain
Material As we consider the material as a continuous medium, we use Representative Volume Element (RVE) [17] method to get the average property (〈߯〉 బ ) of a multiphase system which can be calculated by: where ߯ ఈ is the intrinsic property of phase ߙ , ݃ బ ఈ the volume fraction of phase ߙ in domain ‫ܦ‬ and ߙ could be any phase in Table 2. In the following, the notation ݃ will refer to the volume fraction of ‫ܤ‬ in ‫.ܣ‬

Mass conservation
Mass conservation should be respected in both the possibly multiphasic material and gas domains: where ߩ is the local density, ‫ܝ‬ the velocity deduced from NS equation. Multiplying the first equation by 1 − ℋ and the second by ℋ and using mass conservation along the interface [17], we get the global mass conservation: In the classical velocity-pressure formulation, we express the divergence of velocity and couple it with the momentum equation. In fact, from Eq.(5), we can deduce: with ߜ = ௗℋ ௗట is the Dirac function and n is the unit vector normal to the interface. Here 〈ߩ〉 భ is considered to be time and space independent. The velocity ‫〉ܝ〈‬ బ and ‫〉ܝ〈‬ భ may be assumed to be equal to ‫ܝ‬ in the transition zone with ߰ ∈ [−߳, ߳].

Energy conservation
The energy conservation in the whole system can be expressed by: where ܶ is the temperature, ℎ the specific enthalpy (J/kg), ߣ the thermal conductivity and ‫̇ݍ‬ the volume heat source. This is a non-linear equation as all properties including ߩ, ߣ, ℎ may depend on temperature. Therefore these parameters are tabulated in an input file of the solver [15]. As alumina and zirconia ceramics are transparent to the laser, dopants such as carbon particles are added in powder in order to augment heat absorption [18]. Consequently, the absorption varies in different areas. For example, the absorption of substrate is much lower than that of powder as the former is not doped. The absorption of consolidated material (melted from powder) depends on the atmosphere. It is low in air as dopants are combusted. The interaction between laser and particles (dopants and powder) includes conduction, convection and radiation, which is very complex to model. However, we aim to get a global heat source model taking into account the local absorption of material. A model is developed based on the Beer-Lambert law [18] [19], assuming: ‫,ݔ(ߔ‬ ‫,ݕ‬ ‫ݖ‬ + Δ‫)ݖ‬ = ‫,ݔ(ߔ‬ ‫,ݕ‬ ‫݁)ݖ‬ ିఈ௭ ( 9 ) where ߔ is the flux (W/m ଶ ) which penetrates the material, ‫ݖ‬ the penetration depth in the laser beam direction (Figure 2) and ߙ the local absorption coefficient. On the other hand, we can get the volume heat source by calculating the attenuation of flux in the propagation direction: Inserting Eq.(9) into Eq.(10) and considering the first order development of ݁ ିఈ௭ , we obtain a differential equation: − ݀ߔ ‫ݖ݀‬ = ߙߔ ( 11 ) with boundary condition: which is supposed to be a Gaussian distribution expressed in the local reference of the laser. Here ܴ denotes the reflection coefficient at the powder bed surface, ܲ the laser power and ܴ the beam radius. Solving Eq. (12), the heat source model can be deduced: Eq. (13) shows that this heat source model takes into account the attenuation of flux (1 and 2 in Figure 2) in the laser beam direction (uniaxial integration of α in exponent) and the local absorption of material (coefficient α).

Momentum conservation
In the considered SLM process, there exists 50% porosity in the powder bed at initial state. When the powder is heated and melted, the consolidation due to melting and the evacuation of gas will result in significant shrinkage and the displacement of gas-material interface represented by the LS with ߰ = 0. Therefore, the momentum conservation should take into account compressibility: where ‫ܝ‬ is the velocity, the volume force which includes gravity and surface tension, and the stress tensor which respects the compressible Newtonian constitutive law: ‫ݎݐ‬ ቀࣕ̇ቁ = ∇ ⋅ ‫ܝ‬ = ߠ̇ ( 16 ) where ‫‬ is the pressure, ߟ the dynamic viscosity, the deviatoric stress tensor of , the identity tensor and ࣕ̇ the strain rate tensor. Another important driving force of the interface displacement results from surface tension (N/m ଶ ): ௦ = ‫ܖߢߛ‬ ( 17 ) where ߛ is the surface tension (N/m), ߢ = −∇ ⋅ ‫ܖ‬ the average curvature and ‫ܖ‬ the unit vector normal to the interface. In the present model, the tangential force resulting from the variation of ߛ along the interface (Marangoni effect) is not taken into account.
Numerical methods differ by the formulation of ‫ܖߢ‬ [20]. The explicit method takes it at the previous time step while the implicit method takes it at the present time step, providing more stability. In fact, ‫ܖߢ‬ can be related to the interface position ‫܆‬ by a surface Nabla operator ∇ ௦ : ‫ܖߢ‬ = ∇ ௦ ⋅ ∇ ௦ ‫܆‬ ( 18 ) where ∇ ௦ ‫܆‬ = ‫܆∇‬ − ቀ∇‫܆‬ ⋅ ‫ܖ‬ቁ ⊗ ‫.ܖ‬ We can express ‫܆‬ by the forward Euler method:  ( 20 ) It can be seen that compared with the explicit one, implicit method has an additional term related to the surface diffusion of velocity, which introduces additional stability of resolution. In the level set context, the surface force ௦ must be multiplied by the discrete Dirac function to get a volume force contribution , according to the CSF method proposed by Brackbill et al. [21]. Following Hamide [22], such an implicit formulation was implemented by Khalloufi [23] in the numerical package CimLib used in the current study.
In the present work, we consider that surface tension appears when liquid and gas are present at the same time. This means that it is applied not only at the gas-liquid interface but also at the powder-liquid interface. The former is just one part of the interface ߰ = 0, while the latter interface is immersed in domain ‫ܦ‬ so we will detect it by the contour ݃ బ = 0.5 where gas is present. Note also that the evacuation of gas from the powder bed may affect the flow due to the shear stress (hydrodynamic drag) and possible additional normal load. These effects are not considered here, together with possible material vaporization.

Material properties
In order to calculate the average property 〈߯〉 బ of domain ‫ܦ‬ , we should know the intrinsic property ߯ ఈ of each phase ߙ and its volume fraction ݃ బ ఈ which may be a function of temperature and chemical composition.

Fraction evolution of zones, structures and phases
The melting temperature is constant (1860°C) for the eutectic composition of Al ଶ O ଷ − ZrO ଶ system. It results in a discrete change of enthalpy with respect to temperature, which should be avoided from the numerical point of view. In the present work, a large melting interval (1800 °C to 1900 °C) is taken as shown in Figure 6. Correspondingly, the transformation from powder (ܼ ) to compact medium (ܼ ଵ ) begins at the same temperature (1800 °C) as material is heated rapidly and the consolidation due to sintering is negligible. However, as powder has a low conductivity, the thermal gradient in powder around the laser is very significant. It means that a very fine mesh should be used in this area in order to well capture the shrinkage.
Another choice is to assume a greater temperature interval for ܼ to ܼ ଵ transformation (Figure 7), which is equivalent to suppose a sintering step before melting. This transformation is non-reversible, so arrows are used to indicate heating and cooling path. In powder, ݃ బ ௌ బ = ݃ బ ௌ భ = 0.5 as there is 50% porosity. The volume fraction of each phase in ܵ and ܵ ଶ can be easily deduced according to the eutectic composition and the density of each phase.

Properties
The specific enthalpy of Al ଶ O ଷ and ZrO ଶ phases as function of temperature is calculated by the model in [24]. Specific enthalpy out of the valuable interval is extrapolated. The specific enthalpy of liquid phase is averaged between liquid Al ଶ O ଷ and liquid ZrO ଶ by their mass fraction. The specific enthalpy evolution of material (Figure 8) is tabulated for the thermal resolution (Eq. (8)). The absorption of material domain ‫ܦ(‬ ) is averaged between powder (ܼ ) and compact medium (ܼ ଵ ) by: where 〈ߙ〉 బ = 70 mm ିଵ and 〈ߙ〉 భ = 60 mm ିଵ are supposed [18]. The gas domain ‫ܦ(‬ ଵ ) is considered as transparent to laser. The coefficient of surface tension is 0.55 N/m for gas-liquid interface. The density of each phase is given in Table 3. Table 3. Density of each phase [25] ݉ ‫ݐ‬ ܽ ݈ ݃ ߩ (kg/m ଷ ) 5680 5680 3800 4405 1.3 Thermal conductivity is a crucial parameter to the temperature distribution in SLM. The conductivity of powder is much lower than that of compact medium. As a consequence, most part of the heat is extracted by the substrate. Similar to the absorption 〈ߙ〉 బ in Eq.(21), we average the conductivity of powder (ܼ ) and compact medium (ܼ ଵ ) to get 〈ߣ〉 బ and 〈ߣ〉 భ : The conductivity of Al ଶ O ଷ could be found in [26] with expression ൫ܶ ∈ [25, 1300] °C൯: It is extrapolated out of the temperature interval. The conductivity of ZrO ଶ (2 W/(m K) [27]) is considered constant. Several models exist for the calculation of powder conductivity and we have taken the Zehner-Schlunder model [28], which considers the conductivity of gas (ܵ ଵ ) to be dominant compared with that of solid (ܵ ). The conductivity evolution of powder is shown in Figure 9. Figure 10 shows the conductivity of material domain ‫ܦ(‬ ), which follows the red line (same as 〈ߣ〉 బ when there is no ܼ ଵ ) when the material is heated. Once it is melted and then cooled, it follows the blue line (〈ߣ〉 భ ).   The dynamic viscosity of liquid Al ଶ O ଷ and ZrO ଶ can be found in [29]. The dynamic viscosity of solid is taken as 1000 Pa ⋅ s. For compact medium (ܼ ଵ ), 〈ߟ〉 భ is averaged by the geometric relation: The dynamic viscosity of material ‫ܦ(‬ , Figure 11) is then calculated by:

Simulation configuration
As illustrated in Figure 12, the whole system simulated is 0.5 × 0.2 × 0.3 mm ଷ and it is separated by the LS interface at z = 0.2 mm at initial state. Gas is above the interface and material is below. A substrate is located at the bottom ‫ݖ(‬ = 0 ∼ 0.17 mm) and a layer (  The initial temperature is 20 °C. Convection and radiation conditions are applied at the bottom of the substrate with heat transfer coefficient ℎ ் = 1000 W/(m ଶ K) and emissivity ߳ ் = 0.4. Given the small thickness of the substrate in the system, an arbitrary augmented value for ℎ ் is chosen, in order to simulate heat extraction from a real size substrate. Radiation of gas-material interface with ߳ ் = 0.6 is also taken into account. By means of the CSF method, following Desmaison et al. [12], a contribution to the volumetric heat source term ‫̇ݍ‬ in Eq.(8) can be expressed. Other surfaces are adiabatic. The bottom is fixed ‫ܝ(‬ = 0) and the top is a free surface through which gas can enter into the system. The normal velocity along the four lateral faces is supposed to be zero. Blue line (Fig. 10) 27 27 As the powder is almost adiabatic, the conductivity of compact medium (〈ߣ〉 భ ) is very important to temperature field. Other important parameters are the pair of process parameters (ܲ , ‫ݒ‬ ). Therefore 3 tests are compared to investigate qualitatively the influence of ‫ݒ‬ and 〈ߣ〉 భ . Parameters are listed in Table 4.
In order to simulate the effect of possible Marangoni flow on the homogenization of temperature in melt pool, 〈ߣ〉 భ is multiplied by an arbitrary factor 5 when ܶ > 1900 °C. This modification is directly inspired by the modeling of fusion zone in fusion welding [30]. The time step is 1 μs for all tests.

Results and discussion
At first view, the bead becomes narrower when 〈ߣ〉 భ and ‫ݒ‬ increases (Figure 13). Beside the bead, a tough edge is formed and it is still stuck to the powder bed. This is a zone partially melted and then solidified. The dynamic of bead formation is different for the 3 tests. In Test 1, as the conductivity of compact medium ܼ ଵ is very small, the melt pool (T > 1900 °C) expands laterally and is larger than the laser beam diameter. The powder in front and along lateral sides of melt pool is wetted by it. When the compact medium is more conductive, like in Test 2, the bead becomes narrower. The bead formation is then found periodical, as shown in Figure 14. At first (a), the powder is not wetted by the melt pool while it is affected by the laser beam. Latter (b) the affected powder is melted and droplets are formed and spheroidized by surface tension. Then (c), they fall down and join the melt pool. The local volume increase leads to move forward the melt pool, causing an augmented wetting of the powder. This phenomenon results in slight waves on the surface of the melt pool. At last (d), the melt pool comes back to the form in (a) due to surface tension and the waves disappear after solidification. These steps are continuously repeated while the laser moves forward. As seen from Figure 14, the hydrodynamics of the melt pool shows spatial variations with velocity varying between v and 2.5 v (when droplets joint into the melt pool). In Test 3 (Figure 15), the bead is developed continuously and the melt pool appears to be more stable, being always in contact with the powder bed, and not showing a marked periodic regime like in Test 2. However, from time to time, it is disturbed by the formation of droplets. It may be thought that such perturbation is related to the asymmetry and the size of mesh, but this requires further investigation. The temperature distribution and the form of melt pool are compared in Figure 16. In Test 1, heat can't be extracted by the substrate and it is blocked in the melt pool due to the low conductivity of compact medium Z ଵ . The temperature in the melt pool exceeds 4000 °C. The melt pool is very long, large and deep. With the same scanning velocity, Test 2 shows a smaller melt pool with ܶ < 2200 °C as the compact medium ܼ ଵ is much more conductive. Comparing Test 2 and 3, we can see that the melt pool becomes longer but shallower and narrower when the scanning velocity increases. In all tests, the melt pool is lagged behind the laser beam. Thermal gradients are higher in powder than in compact medium due to different conductivities.

Conclusions
In the present work, the 3D finite element modeling of bead formation during SLM additive manufacturing has been presented. The Level Set method is used to track the gas-material interface and the multiphase homogenization takes into account the different phase transformations of the material. As dopants are added to monitor the laser absorption, the heat source is considered on volume rather than restricted to the material surface. Thereby, a heat source model based on Beer-Lambert law is deduced, which takes the absorption into consideration. The thermal resolution is coupled with thermodynamics. The shrinkage during consolidation from powder to compact medium due to the evacuation of gas is modeled by means of a compressible Newtonian constitutive law. Semi-implicit formulation of surface tension is used, which permits a stable numerical resolution to capture the shape of melt pool.
Three tests with different conductivity and scanning velocity are compared to investigate qualitatively their influence on temperature distribution, dimension of melt pool and dynamics of bead formation. Although many parameters (thermal conductivity, absorption, reflection, etc.) need to be confirmed by experimental measurement or inverse method, our model shows an excellent sensitivity to different parameters. Future work will be focused on the calibration of parameters and their influence on multipass deposition.