*2.1. Constitutive Equations and Its Implications on the Modeling*

### 2.1.1. Absorption Coefficient of the Metal Powder

When the laser beam hits the metallic powder, part of the electromagnetic radiation is absorbed leading to an elevation of the temperature in the interaction area. During the time in which the powder is heated but not having experienced melting yet, that part of the radiation that is reflected from a given particle may impact on one of the surrounding ones, which, in turn, will behave in the same way absorbing part of the radiation and reflecting the remaining. The repetition of this process among adjacent particles (see Figure 1) leads to a higher overall absorption coefficient, especially, if compared with the situation of a laser beam irradiating a flat surface where the reflected radiation escapes from the material.

The total absorption, *AT*, can be characterized by the addition of the radiation absorbed by each particle *Ai* plus the effect of the reflected radiation on each particle, when it impacts on the surrounding ones, as the product of the reflection coefficient of the first one *Ri* by the respective absorption coefficient corresponding to the particle which receives this reflection *Ai+1*. According to the Fresnel Equations from the electromagnetic theory [14], the absorption coefficient of radiation when it impacts on the surface of a conducive media depends on the kind of polarization of the radiation, its wavelength, the optical properties of the material (complex refraction index) and the angle of incidence between the propagation direction of the radiation and the normal vector of the irradiated surface, *θ*. While most of the indicated variables can be fixed for a given process configuration, and the complex refraction index of the material for the process wavelength can be taken from [15], the incidence angle of the laser radiation on each particle is a random variable, whose value can vary between 0 and π/2 rad. Since there is no reason to consider any value of the defined range as the most likely, this uncertainty can be solved by taking the average absorption coefficient from all the values of the range, *A*(*θ*<sup>0</sup> *<sup>π</sup>*/2). In this way *Ai+1 = Ai = A*(*θ*<sup>0</sup> *<sup>π</sup>*/2). Additionally, not all radiation reflected from a given particle is trapped by the surrounding ones, and it can be released definitely from the bulk of powder. This phenomenon can be described by a coefficient to retain radiation within the bulk of powder, *η*, which can be related with the thickness of it, *h*. Like this, the multiple reflection process can be described mathematically by a series, as Equation (1):

$$A\_T = A(\theta^0\_{\pi/2}) + \eta(h)R(\theta^0\_{\pi/2})A(\theta^0\_{\pi/2}) + \eta(h)^2R^2(\theta^0\_{\pi/2})A(\theta^0\_{\pi/2}) + \eta(h)^3R^3(\theta^0\_{\pi/2})A(\theta^0\_{\pi/2}) + \dots \tag{1}$$

Equation (1) corresponds to a convergent series that can be expressed by a general term as Equation (2):

$$A\_T = \frac{A(\theta^0\_{\pi/2})}{1 - \eta(h)R(\theta^0\_{\pi/2})} \tag{2}$$

**Figure 1.** Illustration of the way in which the radiation reflected from a given particle may impact on the surrounding ones, in a multiple reflection process, resulting in a high absorption of the radiation supplied by the laser beam.

#### 2.1.2. Thermo-Fluidic Coupling

Once the energy from the laser source has been absorbed by the powder, its temperature quickly increases. At the beginning, the short contact surface of the particles limits thermal diffusivity and the melting of the powder occurs just in the proximities of the laser-powder interaction zone. Once the liquid has been formed, the heat diffusion and wettability conditions improve, and the melt pool spreads beyond the limits of the interaction area [16]. The first phase change that takes places happens, therefore, for that part of the powder that becomes liquid, and the subsequent phase change will happen with the solidification of that liquid to form a consolidated solid. Both transitions, powder to liquid and liquid to solid, are driven by the thermal cycle imposed by the scanning of the laser source (fast heating and fast cooling respectively), although, each of them takes places under specific circumstances. The material which experiences melting during the heating stage evolves from being disaggregated particles with very low cohesion force among them and no continuous surface (each particle has its own one), to a continuous medium where exists friction among molecules characterized by the dynamic viscosity of the liquid metal, and, at the surface, cohesive forces associated to the surface tension of the liquid. The shape of the melt pool is configured, therefore, by the combined effect of the own weight of the melt pool (hydrostatic pressure), the influence of dynamic viscosity and surface tension. This combination of loads leads, as will be seen with the experimental and numerical results, to a geometry of the melt pool tending to acquire drop shape.

The solidification of the liquid, during the cooling stage, takes place by the increase of the cohesive forces among molecules of the material during the consolidation phenomena. In this way, the shape acquired by the melt pool tends to be maintained despite the disappearing of the effect of the balance between the surface tension and the hydrostatic pressure.

Figure 2 schematically represents the traversing of the laser beam across a given section for different temperatures, T, in order to illustrate the preceding ideas. The final height of the consolidated material, as well as its drop shape are the consequence of several physics phenomena. In the first place, the poor thermal conductivity of the powder makes the energy supplied by the laser beam concentrate intensely around the laser-powder interaction area. It makes that in that region the melting point is exceeded, not only by the powder and substrate which are strictly under the laser beam, but also by the powder around it. It makes that, in the end, the volume submitted to melting is larger than just the volume of powder strictly under the laser beam. Once the volume of heated powder becomes liquid, its surface tension forces it to evolve form a volume with the shape of a rectangle in the powder be to a taller drop of liquid metal, whose height is dependent mainly of the time that the material evolves in liquid state. Additionally, given the low amount of energy needed to melt isolated particles, those of them around the melt pool are wetted by it and consequently melted, being incorporated into the melt pool, contributing significantly to the amount of liquid material. The concentration of energy under the laser-powder interaction area must also allow for melting the substrate to ensure a proper adhesion of the consolidated material to it. If the melting of the substrate does not happen, the melt pool does not penetrate it, and, although a consolidated ribbon might be formed, it is not welded to the substrate. During the melting of the powder, the air contained initially in the pours is released, making the density of the domains increase. Again, during the cooling of the material some effect of shrinkage happens, which is also reflected in the evolution of density. This effect during cooling is also considered in Figure 2.

**Figure 2.** Illustrative representation of the consolidation of the material during the traversing of a given cross section by the laser beam.

The thermal part of the calculations is governed by the heat conduction equation in a fluidic media, Equation (3);

$$
\rho \mathbb{C}\_p \frac{\partial T}{\partial t} + \rho \mathbb{C}\_p \mathbb{u} \cdot \nabla T = \nabla(\kappa \nabla T) + Q. \tag{3}
$$

In Equation (3) the term *Q* represents any volumetric energy exchange that takes places in the studied bulk material. In this context, the latent heat of the phase changes and the metallurgical transformations are the inputs for that term. Equation (4) introduces the laser source as a surface boundary condition:

$$-\mathfrak{n} \cdot (-\mathfrak{x} \,\nabla T) = Q\_b \tag{4}$$

In Equation (4) *Qb* represents the irradiance of the laser beam.

The fluidic evolution of the melt pool, when the material is in liquid state, is governed by the Navier–Stokes equation, as Equation (5) indicates:

$$
\rho \frac{\partial \mathbf{u}}{\partial t} = \nabla [-pI + \mu \,\nabla \mathbf{u} + \mu (\,\nabla \mathbf{u})^T] + \rho \,\mathbf{g} \,\mathbf{z} \tag{5}
$$

In Equation (5), the last term, *ρgz*, includes the weight of the liquid as a load in the balance of forces. *ρ* is the density of the material, *g* is the gravity constant and *z* is the height under the surface of the fluid. The tangential force at the surface of the fluid *Ft* because of the surface tension *σ* is described by Equation (6):

$$-\nabla \sigma + \sigma (\nabla \cdot \mathfrak{n}) - p\_{\text{ext}} \mathfrak{n} = \mathfrak{n} \cdot \mathbf{F}\_{\mathfrak{t}}.\tag{6}$$

In Equation (6) the term *pext* includes the effect of the external pressure (the atmospheric pressure at room temperature) on the surface of the fluid.

Figure 3 shows a flux diagram of the thermo-fluidic coupling. It happens, therefore, due to the phase and temperature dependence of some critical properties of the matter, such as dynamic viscosity, *μ*, surface tension, *σ*, and density, *ρ*. The specific capability, *Cp*, and thermal conductivity, *κ*, which are associated exclusively with the thermal balance, also show phase and temperature dependence. For numerical purposes, a variable to identify the phase of the matter is defined, Ψ (Phase Variable), its value is set to 1 if the material has exceeded the melting point, and 0, if it has never happened up to the analyzed instant.

**Figure 3.** Flux diagram of the constitutive equations, boundary conditions and initial conditions of the thermo-fluidic coupling.

Figure 3 shows the coupling between the thermal calculations associated to the heat conduction equation and the laser beam irradiance *Qb* as a surface boundary condition (block to the left) with the fluidic calculations by means of the Navier–Stokes equations submitted to surface tension, σ, and the external pressure, *Pext*, at the boundaries of the domain (block to the right) [17]. The top block in Figure 3 contains the initial conditions with the process parameters included. Some of them, such as the power of the laser *P*, the beam diameter *ø* (a Gaussian beam is considered), the scanning speed of the laser *V*, the heat capability *Cp* and thermal conductivity *κ* are associated with the thermal calculations. Other properties, such as dynamic viscosity *μ* and surface tension *σ*, participate only in the fluidic

calculations. In the case of density of the material, *ρ*, it influences the thermal balance since thermal diffusivity *χ* is obtained by Equation (7):

$$\chi = \frac{\kappa}{\rho \mathbb{C}\_p} \tag{7}$$

and density also appears in the Navier–Stokes Equation. The velocity field of the domain, *u*, represents the three spatial components of the speed of any point of the domain. At the beginning it is set to zero, since there is no movement of the particles in the bulk of powder. On the one hand, the motion of the fluid affects the thermal balance; on the other hand, it is the main output of the fluidic calculations. It allows for determining the movement of the fluid, and therefore, the evolution of its geometry. At the first step of calculation, the properties of the material derive from the initial temperature, *T0*, and initial phase, Ψ *= 0*, since the material has not experienced melting yet. The effective thermo-fluidic coupling is driven, as is implicit in Figure 3, by the thermal cycle, and therefore, by the heat calculations. Once the initial step has been taken, the rising of the temperature determines the turning of the affected points as Ψ *= 1* (the material has molten in these points). The temperature and phase dependence of dynamic viscosity and surface tension regulate, in turn, the flowing of it, associated to the heating, and, at the same time, the evolution of density, from the low apparent density of the powder to the higher density of the fluid, is able to represent the releasing of the interstitial air of the powder, and the appearance of the liquid as a continuum medium.

During the cooling stage, the use of phase-dependent thermo-fluidic properties, along with the phase variable, Ψ, determines if the material must cool, either, as unmolten powder or as a rigid solid.

The theoretical basis for the consolidation of the material has been established: the driving force to turn the powder into liquid is the thermal cycle associated with thermo-fluidic properties dependent on the temperature and the phase of the matter. This dependence also determines the phenomena during cooling, giving rise to a rigid solid in those regions where the melting point has been exceeded. In this context two, challenges arise for the modelling work. From a physics point of view, the behavior of the thermo-fluidic properties must be adjusted considering the requirements of temperature and phase dependence. From the numerical point of view, the distinction between molten and unmolten region by means of a discrete variable, Ψ, has to be implemented on the frame of a moving mesh, capable of representing the movement of the melt pool.

#### 2.1.3. Numerical Tools to Carry Out the Thermo-Fluidic Coupling

The main goal of the thermo-fluidic approach is the calculation of the shape acquired by the liquid material under the effect of the surface tension, gradient of the surface tension, the hydrostatic pressure, the variations of dynamic viscosity and density, as well as maintaining that complex shape during the solidification of the material. It has some important implications for the numerical calculations.

Figure 2 in the previous section illustrates the cross section of a drop formation during the scanning of the laser beam. It can be seen that, starting from a flat contour representing the powder bed, the domain has to evolve to a convex contour several times higher than the original thickness of the powder bed.

The Arbitrary Lagrangean–Eulerian formulation of the Finite Element Method [11] is able to deal with the movement of the mesh and can be used to adapt it to the deformation imposed by the fluidic loads. The basis of this numerical methodology lies in the use of different coordinate systems to have a reference of the movement of the mesh. While the Eulerian part of the calculation uses a fixed reference system, the Lagrangean part uses a reference system associated with the points of the material which moves with them.

Nevertheless, the motion capability of the nodes of the mesh is not enough to represent high complex geometries. The movement of the nodes of a given element with different directions may lead to it having a high level of distortion.

According to the general theory of the Finite Element Method [18], a high degree of distortion of the elements negatively affects the accuracy of the calculations and may lead to the impossibility of performing the convergence of the solution. It happens due to the fact that the values of the analyzed function are calculated, in the interior of each element, by interpolating the values of the nodes. If the aspect ratio of the element is too high, the result of that interpolation is wrong.

The use of the Automatic Remeshing technique allows the problem of the distortion of the elements to be managed. It is a feature directly available in the software used for the simulations; COMSOL MultiphysicsTM. In this way, the distortion is calculated internally by the program as a positive dimensionless number for each element. The value of this number for an undistorted element is 0. It can rise and, even descend, following the geometrical evolution of the element. When the distortion of any element reaches a critical value, the calculations are stopped, and a new mesh is built on the basis of the deformed geometry at the current simulation step. Figure 4 illustrates the simulation of the formation of a drop as a consequence of the melting of the powder. The initial mesh (left) resembles the flatness of the bed of powder. From that moment the simulation proceeds and the droplet appears progressively (center). At a given instant the distortion reaches the maximum permitted value. In that instant the numerical calculation stops, and a new mesh is automatically built by the program (right), if not exactly equal to the deformed geometry, optimized according to the algorithms of the program, and the initial user-defined mesh size.

**Figure 4.** Process of an Automatic Remeshing. Departing from the initial mesh (**left**), the drop formation starts (**center**), and, when a high value of distortion in any element is reached, a new mesh is built, based on the deformed geometry (**right**).

The maximum value of the distortion at which a new mesh has to be done is user-defined. The procedure to adjust this value starts by attempting to carry out the simulation without limiting the distortion of the elements. The simulation finds problems in making the solution converge with a value of distortion about 20, and the calculations can no longer carry on. From this evidence it may be thought that by limiting the distortion to a very low value the convergence of the solution is guaranteed. Nevertheless, low values of distortion, such as 1–2, may be quickly reached and forced to stop and resume the simulation very frequently, with the associated wastes of time and memory. Considering the aforementioned restrictions, a maximum value of 5 has been used as a criterion to proceed with the automatic remeshing.

### *2.2. Thermo-Physical Properties for the Thermo-Fluidic Coupling*

Thermal conductivity shows a high phase-dependent behavior. [19] highlights the influence of the limited contact surface among the particles or the powder bed as a limiting factor for the heat to diffuse. In the work referred to, where the AISI 316L steel is considered, the thermal conductivity of the powder at room temperature is significantly shorter than the conductivity of the same material in solid state. The present work proposes the use of different temperature functions depending on the absolute maximum temperature reached for each region of the domain. If a given region has reached, at any time, a temperature corresponding to a melting condition, its phase variable is set to, Ψ = 1. It determines that its cooling happens by a function different from the one during the heating, which, instead of returning to the initial conductivity of the powder, evolves along the conductivity corresponding to the cooling of the solid, considering, even, the different solid state metallurgical phases, if it is the case. Figure 5 shows the temperature and phase dependence of thermal conductivity.

**Figure 5.** Temperature and phase dependence of thermal conductivity.

An important feature of the proposed behavior is the transition between either, powder to liquid or liquid to solid. It considers conditions far from the thermodynamic equilibrium associated to the high speed of the heating induced by the laser in the material, and the subsequent cooling, which is something considered in some references for the case of high power laser treatments [20]. In this way, the change of state, instead of being quasi-instantaneous, as for a situation of equilibrium, takes place along a temperature interval following a sigmoid behavior. It is reflected on the thermo-physical properties.

Density, considering the same hypotheses as for thermal conductivity, shows an intuitive behavior. It evolves from the apparent density corresponding to the powder of metal, to the density corresponding to the metal liquid. During the cooling, those regions which have molten, use the evolution of density with temperature corresponding to a solid. Figure 6 illustrates the phase and temperature dependence of density.

**Figure 6.** Temperature and phase dependence of density.

Dynamic viscosity plays a crucial role in the evolution of the shape of the material treated. The value of the dynamic viscosity of a liquid steel can be found in [16]. It can be maintained constant along the liquid state. Dynamic viscosity can also be used to emulate the behavior of a solid during the cooling of the material. By associating the dynamic viscosity at room temperature with a value high enough, the shape acquired by the melt pool during the time that the material is in liquid state, can be maintained when it cools. This methodology is proposed by [21] for the keyhole welding process, where a value of 100 Pa·s is considered to be enough to reproduce the behavior of the solid metal. The value of dynamic viscosity to represent the behavior of the powder of metal is the variable surrounded

by a higher level of uncertainty. In this work it has been estimated by considering the size and the dimensions of the droplet formed. The formation of it takes place in a volume where the droplet is surrounded by powder. Several experiments have been carried out by using a powder bed with a relatively high thickness of 1 mm. The objective of using such a large layer lies in investigating whether a big volume of powder tends to limit somehow the capability of the melt pool to acquire a free shape. If this is not the case, it can be concluded that powder presents no significant wear for the development of the melt pool, and, therefore, its dynamic viscosity is negligible. The fundamentals of this analysis lie in the comparison of it with the re-melting test where a laser beam is applied on a solid metallic plate [22]. When this is done, despite the material reaching the liquid state, no droplet formation happens. It can be associated with the limitations that the solid metal surrounding the melt pool imposes on the free movement of it.

Figure 7 contains the image of the cross section of several ribbons obtained after the scanning of the laser beam on a powder layer of 1 mm. In all the cases the power of the laser beam was 4000 W and the scanning speed of 400 mm/s, 600 mm/s and 800 mm/s.

**Figure 7.** Cross sections of the consolidated material made with a layer of 1 mm of powder with 4000 W, 400 mm/s (**left**), 600 mm/s (**center**) and 800 mm/s (**right**).

From the results in Figure 8 it can be assumed that the powder around the melt pool, despite having a high thickness, it has not prevented the liquid material from adopting the corresponding drop shape. Consequently, the dynamic viscosity of the metallic powder is considered negligible in comparison with the value of the liquid. Figure 8 contains the proposed phase and temperature functions for dynamic viscosity.

**Figure 8.** Phase and temperature dependency of dynamic viscosity.

Finally, the surface tension is the last variable that is introduced as critical for the definition of the shape of the melt pool. In the liquid state it shows dependence with temperature as line function with negative slope. The particular value for the working material is reported by [8]. The consideration of this dependence with temperature in the liquid state plays a crucial role in the definition of the melt pool shape, since it is related to the so-called Marangoni Convection. Under the influence of this phenomenon, the liquid metal flows from the regions with low surface tension to the regions where the surfaced tension is higher. In the case of the metal powder, since there is no continuous surface, it can be considered that the surface tension shows a negligible value. For the case of the solid metal the attribution of a value for the surface tension makes no sense. The material is going to be frozen under the effect of a high value of dynamic viscosity and no tangential effects on its surfaces can affect the shape of it.

Conclusively, the important aspects of the surface tension are the transition of the material from powder to liquid, and the temperature dependence of it in the liquid state. Figure 9 shows the phase and temperature evolution of the surface tension. Reference [23] highlights the influence of the surface tension in the Selective Laser Melting Process by the effect of the capillarity force on the melt pool.

**Figure 9.** Phase and temperature dependency of the surface tension.
