Evolution of Melt Pool and Porosity During Laser Powder Bed Fusion of Ti6Al4V Alloy: Numerical Modelling and Experimental Validation

Laser powder bed fusion (L-PBF) is one of the most promising additive manufacturing technologies for metals. In this work, the discrete element method (DEM) was used to reproduce a powder bed of particles distributed in a random way to be as close as possible to reality. Single and multiple scan tracks simulations were performed for Ti6Al4V alloy using a commercial CFD software, FLOW-3D AM®. The output from the numerical simulations was elaborated to obtain shape and size of melt pools, morphology of scan track surfaces and porosity content. In particular, a specific model was used in order to predict air entrainment in the melt pool and, therefore, to estimate gas porosity content, as an innovative approach to predict such defects. Results from simulations were compared with experimental data from Ti6Al4V samples produced by L-PBF in terms of melt pools size and morphology, as well as density. The good agreement between calculated and experimental results indicates that simulation of L-PBF can represent a powerful tool not only for the optimization of process parameters, but also for the prediction of porosity level.


Introduction
Laser powder bed fusion (L-PBF) is an already established additive manufacturing technique, in which a laser selectively melts a pre-deposited bed of microscopic metal powders [1]. The localized laser irradiation results in a micro-sized zone of molten material, i.e. the melt pool. As the laser beam moves on, the melt pool solidifies very rapidly, forming a track of solid material. These steps are repeated until the part is gradually built in three-dimensions, layer by layer. Many industries, such as aerospace or biomedical ones, have interests in employing AM to produce application-tailored components [2].
Despite of the undeniable advantages of this technique, such as the capability to produce complex parts within an acceptable timespan and with low material waste, quality is recognized as a critical issue for AM parts [3]. Variability in component quality comes from the rapid solidification typical of this process [4], as well as defects appearance during metal AM [5]. Figuring out the complex phenomena that occur during metal AM processes is very challenging since real-time monitoring of the process at such mesoscale and temperature is very hard. Furthermore, trial and error method to reach optimum part quality is usually time-consuming and costly.
In this regard, simulations of AM processes have been recently pursued to optimize process parameters and to predict the formation of defects occurring in and close to the melt pool and to find out a way to avoid their formation [6,7]. Thus, process optimization could be achieved without in-line production trials, leading to considerable economic advantages and timesaving. Various studies dealt with the simulation of the formation of defects such as balling and spattering [8][9][10]. Porosity is also one of the most harmful defects occurring in AM components [11] since pores act as stress concentrators so that mechanical properties, in particular strength and fatigue resistance, are decreased [5]. Therefore, many studies have been carried out on this topic. Bayat et al. [12,13] focused on the creation, evolution and disappearance of keyhole and keyhole-induced porosities. Khairallah et al. [14] used a three-dimensional high-fidelity powder-scale model to analyze the influence of recoil pressure and Marangoni convection to define pores formation mechanisms, in particular the effects of depression collapse on porosities and material spattering. Martin et al. [15] explored the collapse of deep keyhole depressions and the consequent trapping of inert shielding gas by means of multi-physics simulation and X-ray imaging. In a study conducted by Yang et al. [16], the presence of spatters was identified as responsible for an increase in pores formation. Khairallah et al. [10] also studied the transition to keyhole induced defects as a function of laser power. Other authors conducted a numerical analysis about keyhole behavior and related pore formation during L-PBF [17]; they found that keyhole became stable at high power thus reducing the risk of pores appearance. Tang et al. [18] focused on lack of fusion porosity simulations; they determined the values of process parameters, such as layer thickness and hatch spacing, at which porosity increases. Similarly, Teng et al. [3] performed a layer-by-layer analysis to predict lack of fusion porosities in CoCr parts. Vastola and co-workers [19] developed a numerical model for porosity prediction at high beam energy regime when transition from conduction mode to keyhole mode occurs.
It appears that several works have focused on the prediction of keyhole or lack-offusion porosities. However, various gas entrapment phenomena usually occur during the evolution of melt pool and the formation of gas porosities is extremely likely to happen and, at the same time, difficult to be avoided. Nevertheless, little attention has been devoted to numerical simulation of porosities due to gas entrapment. In fact, this is more complex than for lack-of-fusion porosities, that can be easily identified by the mesh cells because of their usually larger sizes than gas pores. However, the detrimental effects of the latter have to be properly considered during design and optimization of a new component, particularly for structural applications, where heat treatments are mandatory, but they can lead to an increase of voids content [20,21]. In this regard, numerical simulation appears extremely useful for the prediction and analysis of the formation of gas porosities. Hence, in this paper the numerical simulation was used to study the fluid flow of melt pools and the evolution of scan track profile in a L-PBF process, as well as to predict the level of porosities due to gas entrapment. Numerical simulations were carried out by means of a powder-scale CFD commercial software, FLOW-3D AM®. To validate the results from simulation, experimental analysis was carried out on samples made of Ti6Al4V alloy produced by L-PBF. Model validation was achieved by comparing melt pool size and surface morphology along the melt track with experimental results. In addition, the use of a specific model for air entrainment prediction was proposed to estimate the level of porosities due to gas entrapment.

Samples Production and Characterization
Cubic samples (10 mm 3 ) were produced from Ti6Al4V inert gas atomized powder (Table 1) using a M2 Cusing (Concept Laser, Germany) L-PBF system, with an island scanning pattern. The particle size distribution (PSD) of the used powder is reported in [22]. The parameters used to fabricate the samples were 300 W laser power, 1750 mm/s scanning speed, 50 μm layer thickness, hatch spacing of 75 μm. The process was conducted in a chamber filled with Argon gas to minimize oxygen pick-up to < 0.1%. In this case, no supports were needed to assure the stability of the cubes.
The top surface of the obtained sample was observed by LEICA DMS300 digital microscope and LEO EVO 40 scanning electron microscope (SEM) to characterize the morphology of the scan tracks. Subsequently, the cubic sample was cut to observe the cross section of scan tracks on the top layer. These were considered as representative

3
and directly comparable to the results from numerical simulation of single and double scan tracks in terms of size and morphology of the melt pools. A precision cutting machine was used to reduce any alteration of the surface and of microstructural properties of the sample. After cutting, the sample was mounted in acrylic resin and mirror polished. Chemical etching was performed using Kroll's reagent for 30 s to reveal material microstructure and melt pool boundaries. The sample was observed by means of LEICA DMI5000 M optical microscope. Analysis under polarized light was carried out to identify the grain structure of the material. Width and depth measurements for twelve melt pools were carried out and average values, together with standard deviations, were calculated. Image analysis was carried out on four micrographs at 500x magnification to measure the average size of porosities. Density measurements were carried out according to the Archimedes' method. Weight measurements both in air and liquid (water) were performed using a Gibertini E42-B weight scale. The density was calculated as follows.
Subsequently, the corresponding porosity was calculated using the theoretical density of Ti6Al4V alloy ( 0 = 4.43 g∕cm 3 ).

Numerical Simulation
Numerical simulations were carried out using a powder-scale CFD commercial software, FLOW-3D AM®, to reproduce the L-PBF real process by setting the same parameters used for samples production. First, the discrete element numerical method (DEM) of FLOW-3D AM® software was used to accurately simulate the formation of the powder bed.
Subsequently, the weld module of the CFD simulation software FLOW-3D AM® was used to properly simulate the interaction between laser and powder particles. The meso-scale CFD model is based on the solution of mass, momentum and energy conservation equations (see the Appendix). Furthermore, it includes all relevant physics to capture melt pool dynamics, melting and porosities formation, including equations to properly evaluate the Marangoni effect and the recoil pressure of evaporation [23]. Details on the equations specifically solved are summarized in the Appendix ("Governing equations and thermophysical properties").
In this regard, the innovative aspect of this study is the investigation of the reliability of the simulation in predicting the level of porosity by means of a specific parameter indicated as air entrainment.
Considering the first step of creation of the powder bed, the PSD of the Ti6Al4V alloy powders used for sample production [22] was discretized as requested by the DEM interface. The maximum number of intervals allowed by the DEM module (i.e. nine intervals) was used to describe the PSD (volume percentage as a function of particle diameter). As shown in Fig. 1, this resulted in a reduced interval in terms of particles diameter, ranging from 18 to 61 μm. The volume of particles corresponding to diameters outside this range (lower or larger diameters) was distributed to maintain a good match with the experimental PSD. Consequently, slightly lower values of D 50 and D 90 if compared with experimental results were obtained. Despite these minor differences, a good discretization of the experimental PSD was achieved.
The so-modelled powders were freely dropped on a substrate by means of DEM module. Powders were spread over a substrate of 2.8 × 1.5 × 0.5 mm 3 . A roller was modelled to create a layer with thickness of 50 μm. Powder packing density was calculated by considering particle size distribution and layer thickness.
Powder bed coming from DEM simulations was converted into a .stl file and imported into WELD interface. A solid substrate made of Ti6Al4V was then created. Mesh cell size was set to 5 μm in order to properly discretize the radial gaussian distribution of laser heat flux, as evident in Fig. 2. In addition, this mesh cell size can well intercept the particles in the powder bed since it is significantly lower than the minimum particles diameter (i.e. 18 μm). Therefore, the mesh cell size is adequate to capture the interaction between the powder and the laser, as well as to intercept the volume of material involved in the phase changes from solid to liquid state and vice versa. A finer mesh would result in longer calculation time without providing additional insight into the described phenomena. Table 2 summarizes the process parameters for the welding simulation. Laser beam was set as cylindrical so only the radius is used to describe the geometry of laser beam. This means that the heat flux follows the Gaussian distribution represented in Fig. 2 for each position along the z axis. This indicates that it is assumed that the power of the laser is independent form the distance between the source of the laser and the powder bed. The method for laser-particle interaction is based on the ray tracing method to consider the multiple reflections of the laser rays. Fluid absorption rate coming from literature for the studied alloy is also reported [24]; it describes the fluid capability of absorbing energy from laser, so it deeply affects melt pool morphology. Physical and thermal properties for Ti6Al4V alloy come from the literature [9] and are summarised in the Appendix ("Governing equations and thermophysical properties"). Figure 3 illustrates the geometry for the numerical simulation of a single scan track in the weld interface.
The schematic of the computational domain used in the present study is shown in Fig. 4, where two mesh blocks are visible. An inner mesh is generated for the volume corresponding to Fig. 3. In addition, an outer mesh block is present only to capture the thermal diffusion. For this reason, this outer mesh is slightly coarser than the inner one.
The central domain dimensions are 2 × 1.5 × 0.4 mm 3 and the choice for a regular grid with a uniform size of 5 μm results in approximately three million cells. The time step was set to be in the order of 10 − 8 s. A wider mesh is used for the surrounding volume to allow calculation of heat exchange with cell size of 25 μm in a shorter computational time.
Boundary conditions are set as described in Fig. 5. The letter P indicates that the boundary conditions are atmospheric pressure and temperature. The letter S allows the continuity between one mesh block and the other one. The letter W indicates where there is a solid wall at temperature of 293 K.
Furthermore, the physics-based model is applied to a XY-plane for the simulation of two adjacent scanning tracks (y = 0.0 mm and y = 0.075 mm). The second   scan track is performed with the same scanning direction as the first track (Fig. 3) following a linear scanning strategy. The laser beam is set off after completing the first scan, it moves to go back at the beginning of the track and the second scan is performed right after (uni-directional scanning strategy). Single and double scan tracks have been already proved effective for the study of overhanging situation at high laser scan speed due to melt instability [25]. Furthermore, the melt pool variation and powder bed condition have been investigated using single and double tracks using experiments and simulations to avoid any balling effect [26].
Also, the choice of studying multiple scan tracks with the same orientation is already documented in the literature as a valid tool to study the temperature field development and changes on the melt pool size and depth, which was analyzed alternating three scan tracks on the same orientation [27]. Furthermore, this approach can also be applied to the case of the presence of dual laser beams with linear scan strategy, which was considered to evaluate the reduction of residual stresses distribution [28].
In addition, multi-layer simulations were performed. In this case, the volume obtained after single and double scan track simulations (including unmelted powder and melted and solidified material) was used as substrate. The steps of generation of the volume of particles and the formation of the powder bed for the second layer were repeated following the same procedure described above. The same power, speed and scanning direction were used also for the scanning of this second layer.

Prediction of Porosity Level
As previously reported, an additional parameter is analyzed in the present study to evaluate the level of gas porosity in the material. In fact, a specific in-house numerical model able to consider air entrainment phenomena was used to calculate the volume of gas porosities into the component. The model considers that, when there is a turbulent flow, small liquid elements can be raised above the free surface and they can trap air, that is subsequently carried back into the body of the liquid. This phenomenon depends on the intensity of the turbulences as compared to the effect of stabilizing forces, such as gravity or surface tension. The model is based on a balance between the turbulent kinetic energy per unit volume, that is defined as P t = ρ Q, and the stabilizing kinetic energy per unit volume P d = ρ g n L t + σ/L t .
Where ρ is the density, Q is the specific turbulent kinetic energy, g n is the component of gravity normal to the free surface, L t is the turbulent length scale, σ is the surface tension. Air entrainment occurs when P t is higher than P d , i.e. the turbulent disturbances are large enough to overcome the surface stabilizing forces.
The model defines a scalar variable to calculate the fractional volume of entrained air in the liquid metal. It is assumed that the trapped air is not able to change the volume or the density of the liquid metal. This is reasonable since in the present application the volume of entrained air is sufficiently small. Further details on the model can be found in [29][30][31]. The output is given as entrained air volume and it represents the total volume of air that is trapped in the molten metal because of turbulences in the melt pool. These porosities are usually extremely fine (few microns) and therefore are hardly resolved by numerical methods based on finite elements since they are usually smaller than the mesh cell size itself. Air entrainment model estimates the total amount of gas entrapment into the flow as proportional to the surface area, A s , and the height of the disturbances above the mean surface level, as described in the following equation: Where V is the entrained air volume, t is time, R is the entrainment rate coefficient, A s is the surface area of fluid, ρ is the density of the fluid, g is gravitational acceleration, L t is the turbulent length scale and σ is the surface tension at the liquidus temperature. The entrainment rate coefficient R is a coefficient of proportionality and it must be empirically defined by the user. In the present study, a calibration process of R was carried out in order to define a direct relationship between the entrained air volume V and R for a fixed geometry (corresponding to the single scan track case), material, boundary condition and process parameters. Density and surface tension are from the literature [9]. Sampling volumes were introduced to calculate the percentage of entrained air (i.e. percentage of gas porosity) as the ratio between entrained air volume and the total sampling volume. The entrained air volume was evaluated after complete solidification of the simulated sample.
After this preliminary step, air entrained volume was calculated also for double and multi-layer scan tracks to monitor porosity evolution during subsequent steps of the process. This approach has been already used experimentally to study the porosity formation mechanisms based on the analysis of the layers overlapping with single and continuous tracks [32] and therefore it is considered reasonable also for numerical studies as the present one.
In fact, single and double-scan tracks allowed predicting the level of entrained air of a single layer, while multi-layer simulations enabled to analyze the eventual formation of porosities at the interface between layers.
Since the entrained air volume calculated in this way does not include shrinkage porosities or voids due to incomplete melting of powders, the presence of these macroscopic voids was studied monitoring an additional specific output of the simulation (indicated as volume of voids).

Microstructural Characterization and Density Measurement
The produced cubic sample made of Ti6Al4V alloy is characterized by an island pattern, in which each inner square shows a different laser scanning direction (Fig. 6a). The scan tracks are well visible inside the single square. A detail is shown in Fig. 6b, where the dotted red line indicates where the sample was cut to investigate the cross section.
Considering the microstructural characterization, it is worth reminding that Ti6Al4V is an α + β titanium alloy: Al is an α phase stabilizer while V stabilizes β, leading to an α + β dual phase at room temperature [33]. It has been already proven that the microstructure evolution depends on process parameters, such as layer thickness, laser power and scan speed [18]. As known, titanium alloys microstructure is deeply affected by thermal history and phase transformations depending on cooling rates [34]. In particular, regarding AM products, some studies reported acicular α' martensite in columnar prior β-grains to be a widely diffused microstructure as a consequence of high cooling rates above the martensite start temperature [35]. Low cooling rates can lead to α + β dual phase instead of α', as observed by Simonelli et al. [36].
Observations of the cross section of the cubic sample revealed a lamellar (Widmänstatten) α + β microstructure (Fig. 7a). This is consistent with literature studies as an indication that fast solidification conditions were ensured during the process. Microscope observations under polarized light allowed to reveal columnar grains; they have grown orthogonally to scan tracks because of the thermal gradients occurring during solidification (Fig. 7b), as reported in literature studies [37].
Melt pools belonging to the top layer of the sample were analyzed. Average and standard deviation of width and depth of the experimental melt pools are 84.4 ± 12.6 μm and 68.0 ± 10.0 μm, respectively.
From image analysis, size of porosities in terms of equivalent diameter was measured and an average value of 0.27 ± 0.08 μm was found.
Finally, the average relative density measured on samples resulted 99.26 ± 0.34%.

Numerical Simulations
First, the packing density of the layer of particles was calculated from the DEM module. A packing density of the powder bed of approximately 40% was obtained, showing good agreement with other studies reported in literature [8,[38][39][40]. In Fig. 8, it is possible to observe the evolution of the single scan track with time. The track is quite regular, resulting in a continuous melt pool since balling defects are not present. A reduction in the width is observed during the evolution of the scan track, which accounted for a transition from circular shape to comet-like shape of the melt pool. Laser scanning ensured sufficient penetration depth of the melt volume into the substrate, so that bonding between the powder bed and the substrate was achieved. It resulted in the absence of voids at the interface between powders and the substrate, as shown in Fig. 9. It was also observed that the molten pool surface is depressed underneath laser beam (see red area in Fig. 9a). This is due to the evaporation recoil pressure generating as the material is heated upon its boiling temperature by the intense laser heat input.
Analogously, when simulating a second scan track, a continuous path is formed (Fig. 10). A slight irregularity in the laser track is particularly evident during the second scan track as the area at the maximum temperature (red zone) seems to Fig. 8 Evolution of single scan track profile as the laser moves on. Time is expressed in seconds be wider. This increase in melt pool size can be due to powder bed heating by the laser source. Furthermore, this can also be due to different heat absorption of particles according to their diameter. In fact, laser beam energy follows Gaussian distribution (Fig. 2) so that energy density is maximum at the center of the laser beam and then decreases. Random distribution of particles, combined with a reduction of energy when moving away from the center of the laser beam, can result in a modification of the area at the highest temperature. In fact, this area is enlarged if small particles prevail at the periphery since these are easily and completely melted. In addition, the significant preheating of the first track also contributes to the asymmetric temperature field of the second track, affecting the shifting of the melt region of the second scan track towards the first one.
Subsequently, the two-scans pattern allowed evaluating the accuracy of the overlapping distance, which was shown to be an important parameter for good surface quality and mechanical properties [41]. Yan et al. [42] demonstrated that the hatching distance should be no larger than the width of the remelted region within the substrate rather than the width of the melted region within the powder layer: in this work these conditions are respected, resulting in a proper overlapped zone with no voids or lack of fusion porosities.
Considering melt pool size, these were measured from three different positions of the same scan track. Post processing analysis allowed to directly measure melt pools width and depth (Fig. 11). The obtained average values are reported in Table 3 and compared with the experimental results reported above. It was found that predicted width and depth of the melt pool are in good agreement with experimental results since errors are in the range of 2% and 10% for both depth and width, which are comparable or even lower than those reported in similar studies. For example, Karayagiz et al. [43] found that melt pool width had 15% relative error compared to predicted values, Lee at al. [44] reported a maximum error of approximately 9% for melt pool width, while Heeling et al. [45] observed a mean error of approximately 20% for melt pool size. This is a clear indication of the reliability of the numerical model used.
Furthermore, the aspect ratio of each melt pool, i.e. the ratio between depth and width, was also calculated and its average value resulted to be approximately 0.7. This is consistent with the limited laser penetration and the consequent absence of keyhole mode and of keyhole porosities. Moreover, Tenbrock et al. [46] defined an aspect ratio of 0.8 to be the threshold value above which keyhole mode melting is dominant. Cross sections at different times allow to follow melt pool evolution, up to reaching the classic semi-cylindrical shape after cooling (Fig. 11).
In addition, simulation accuracy was also ensured by comparing surface morphology of scan tracks from simulation and from sample characterization (Fig. 12).
Similarities between simulation and SEM images are evident; V-tracks left by laser are properly predicted by software and their size is in good agreement with experimental values. Moreover, the initial part of scan tracks results wider than  the rest of the track (Fig. 12a): simulations succeed in predicting this variability in an acceptable way since average error is approximately 5%. As above mentioned, two-layer simulations were also carried out to analyze inter-layer porosities. It was found that melt track evolutions of each layer were approximately the same as those described in the single layer simulations. Figure 13 reports temperature distribution for the first scan track of the first deposited layer (a) and the second one (b). As for the single layer simulations, no voids are shown, and a continuous morphology is visible.

Prediction of Porosity Level
As introduced in paragraph Materials and methods -Prediction of porosity level, the air entrainment model requires setting the entrainment rate coefficient R. Therefore, the coefficient R was first calibrated for the simulation of the single scan track in order to understand the relationship between the coefficient and the calculated air entrained volume for a fixed geometry (single scan track), material, boundary condition and process parameters. To this aim, the entrained air content was calculated as a function of various values of entrainment rate coefficient. Based on these results, a linear regression was carried out and the following relationship (R 2 = 0.97) was found: Based on this preliminary calibration, a coefficient value R of 0.12 was applied since it results in an entrained air level of 0.74%, corresponding to the average porosity level measured for the considered samples, under the assumption that no lack of fusion porosities are present, as the worst condition.
Due to the significant difference in terms of the volume considered for the experimental measurements (cubic sample) and the simulation (single scan track) used in this preliminary step, it was necessary to verify the reliability of this value of entrainment rate coefficient for larger volume of material, as in the case of the multi-layer simulation. In fact, as indicated in Eq. (3), the entrained air volume is a function of the surface area. Therefore, double scan tracks and multilayer simulations involve bigger surface areas, and this will affect the calculation of air entrainment volume. Consequently, entrained air volume measurements were carried out also for double scan tracks and multi-layer simulations. A different sampling volume was chosen to properly analyze the material behavior. An average value of trapped air of approximately 0.82% was found for two-scan tracks, while for multi-layer simulations the level of entrained air was predicted as 0.80%.
Since the air entrainment model is not able to predict any macro-voids, a different output was used to identify the eventual presence of lack-of-fusion porosities in order to provide a thorough investigation of the presence of defects in the material. After complete solidification, no macro-voids (i.e. pores with greater size than the mesh cell) were detected in the same sampling volume on the single scan track.
Analogously, the presence of voids between layers was investigated, defining a proper sampling volume in the multi-layer simulation. Again, no macro-voids were detected, confirming the absence of lack of fusion porosities in the interlayer boundary.
Assessed the absence of keyhole and lack-of-fusion porosities, values of air entrained volume obtained from simulations and average porosity from experimental measurements were compared (Table 4). An error of approximately 8% was found in this preliminary study. Despite this error is not negligible, the model appears a promising tool to provide a global index of the content of gas porosities (4) Entrained air % = 2.51R + 0.44 into the component, especially if used as prediction for optimization of process parameters before production.

Conclusions
In this work, L-PBF of Ti6Al4V alloy powder was modelled using FLOW-3D AM® software. The reliability of the numerical simulations was verified by comparing melt pool size provided by numerical modelling with experimental data. An error between 2% and 10% was found for both depth and width, ensuring model accuracy. This was also confirmed by comparing surface morphology of melt tracks provided by simulation with those coming from SEM observations of samples. Then, the entrained air model was used to predict gas porosity level, which provides a global index of voids content allowing to consider even smallest porosities, usually hardly predicted. A first successful estimation of gas porosity content was achieved since discrepancies between simulation values and experimental results, coming from density measurements, was about 8%. Further optimization of this approach is necessary, as development of the present preliminary study, to enhance the effectiveness and reliability of the prediction of gas porosity content for a wider range of process parameters and materials.

DEM Simulation
The discrete element numerical method (DEM) is able to simulate the interaction between particles and between walls and particles in order to reproduce the formation of the powder bed. This approach considers the particles as rigid bodies that do not deform with the shape of perfect spheres, characterized by three velocity components according to a Cartesian system. The walls are considered as rigid surfaces.
This method expresses the contact forces between the particles as a result of two contribution, an elastic force and a viscous dissipation. Furthermore, given the displacement and the relative velocity between the particles, the motion of the particles is simulated based on the calculation of normal and tangential contact forces.
Given two particles i and j, the normal force is calculated as And the tangential force is calculated as Where k is the spring constant, is the displacement, n and t are the viscous drag coefficient, u ij is the velocity vector of each particle, n ij is the distance between the center of the two particles.

CFD Simulations
It is important to remark that FLOW-3D® commercial software for computational fluid dynamic numerical simulation uses area and volume porosity functions based on the FAVOR (for Fractional Area/Volume Obstacle Representation) methods to model complex geometric regions.
In addition, the Volume of fluid method (VOF) [47] is used to track the free surface boundary of the melt pool. In the VOF method, a fluid volume function V F is defined, which is located at the center of the grid. The value ranges between 0 and 1, where V F =0 represents an empty grid (void) and V F = 1 represents a full grid (totally occupied by a fluid). The fluid volume fraction is therefore used to track the interface between two fluids (or between void/gas and the fluid). The boundary conditions are applied to the free surface of the fluid.
The micro-scale CFD model is based on the solution of mass, momentum and energy conservation equations.
The general mass continuity equation for the present study is The equation of motion for the fluid velocity components (u, v, w) in the three coordinate directions are the Navier-Stokes equations, that can be expressed as follows The energy balance can be expressed as  where, h is the enthalpy, t indicates time, � ⃗ v is the velocity of the melt pool, p is the pressure, k is the thermal conductivity, T is the temperature, ̇q is the Gaussian heat source.
The heat source is described with a heat flux function that follows a Gaussian distribution (the energy density gradually attenuates from the center to the outside) along the radial direction as follows: Where Q is the instant surface heat flux, A b is the absorption coefficient, P L is the laser power, R s is the laser radius, and x s and y s are the coordinates of the laser beam center.

Thermophysical Properties
The thermophysical properties of Ti6Al4V alloy used in the present study are summarized in Fig. 14 and Table 5.
Additionally, the data of surface tension of the liquid metal as a function of temperature (reported in Fig. 15) were required in order to consider the Marangoni effect Funding Open access funding provided by Università degli Studi di Brescia within the CRUI-CARE Agreement.
Data Availability Data will be made available on reasonable request.

Declarations
Competing Interests The authors have no relevant financial or non-financial interests to disclose.

Conflict of interest
The authors state that there is no conflict of interest. Fig. 15 Surface tension as a function of temperature for Ti6Al4V alloy