A Temperature-Dependent Heat Source for Simulating Deep Penetration in Selective Laser Melting Process

: Numerical methods for simulating selective laser melting (SLM) have been widely carried out to understand the physical behaviors behind the process. Numerical simulation at the macroscale allows the relationship between input parameters (laser power, scanning speed, powder layer thickness, etc.) and output results (distortion, residual stress, etc.) to be investigated. However, the macroscale thermal models solved by the ﬁnite element method cannot predict the melt pool depth correctly as they ignore the effect of ﬂuid ﬂow in the melting pool, especially in the case of the presence of deep penetration. To remedy this limitation, an easy-implemented temperature-dependent heat source is proposed. This heat source can adjust its parameters during the simulation to compensate for these neglected thermal effects related to the ﬂuid ﬂow and keyhole, and the heat source’s parameters become ﬁxed once the temperatures of the points of interest become stable. Contrary to the conventional heat source model, parameters of the proposed heat source do not require a calibration with experiments for each process parameter. The proposed model is validated by comparing its results with those of the anisotropic thermal conductivity method and experimental measurements.


Introduction
SLM has attracted significant industrial attention due to the possibility of using a large range of materials and producing complex geometrical products [1,2]. However, the SLM process involves the couplings of multi-physics, such as interactions between laser and powder, fluid flow in melting pool, vaporization, and transformation of phase, etc. Experimental measurements allow the establishment of the relationship between input parameters (such as laser power, scanning speed, hatch spacing, etc.) and certain desired properties, while experimental studies are often costly and difficult.
Thanks to the strong development of computing capabilities and numerical algorithms, the numerical method has now become an essential tool to predict mechanical consequences related to the parameters of manufacturing processes. The objective of numerical simulations is to design and optimize manufacturing processes and components; meanwhile, numerical simulation allows a better understanding of the physical phenomena behind the process. According to the modeling scale, the simulations can be composed of two kinds of models: (i) multi-physical simulations; (ii) thermo-mechanical simulations.
The multi-physical models aim to compute more accurately the thermal and velocity fields by couplings of different scales of physics interactions, such as interactions between an electron or laser beam and powders [3], couplings between fluid and solid [4], etc. Numerous high-fidelity multi-physical models [5][6][7][8][9] have been proposed to model the realistic melting surface of the powder layer by tracking the solid-gas interface. Other models [10][11][12][13] present the powder layer as a continuum medium to simplify the complexity, and an effective thermal conductivity and density are defined for the powder layer. The multi-physical simulations at the micro/mesoscale are extremely computationally expensive, and bridging the gap between the grain scale (micrometers) and the part scale (millimeters) is still a big challenge. In spite of the increases in computer efficiency, multi-physical simulation, starting from laser-powder interaction until the final part-scale distortions, is currently impossible.
The current state-of-art shows that thermo-mechanical simulations have been carried out to predict residual stress and distortion in the part scale [14][15][16][17]. As we can see from Figure 1, the thermo-mechanical simulations are performed in two steps. A heat transfer simulation is first carried out to determine the temperature distributions during the process, and the laser-powder interactions are simplified by a 2D heat source or 3D heat source [18]. Then, the temperature distributions are imposed as thermal loading in a mechanical simulation to predict the distortion and residual stresses. However, as the effect of fluid flows in the melting pool and the Marangoni effect are neglected in the thermal model, the prediction of temperature distributions would not be precise enough, which can result in the less accurate prediction of mechanical distortions and residual stresses. Therefore, correctly predicting the thermal results has become one of the key steps in thermo-mechanical simulations. Various numerical methods have been proposed to solve heat transfer equations, such as analytical solutions [19][20][21][22][23] and FE models [24][25][26][27][28][29]. The comprehensive comparison of the analytical and FEM made by Promoppatum et al. [30] shows that good correlations of the melt pool width have been observed between numerical results and experimental measurements, while the comparison of melt pool depth is neglected. Similarly, Romano et al. [27] employed the FE method to simulate the melting pool size in the selective laser melting (SLM) process, and the numerical results tended to overestimate the melt pool width, especially at high energy input. Therefore, a correction of effective absorptivity has been used to be incorporated with numerical results such that the prediction is in better agreement with measurements. Walker et al. [31] have proposed an automated inverse method to calibrate thermal finite element models, which ensure that both the thermal history data and melt pool geometry are predicted with accuracy. According to the author's knowledge, both the analytical method and FE with a 2D heat source can simulate a half-elliptical melt-pool cross-sectional shape. However, in case of deep penetration, the existing heat source in the FEM method must be calibrated with the experimental results for each combination of parameters; therefore, a heat source that does not require calibration with experimental results and is capable of correctly predicting the melt pool's dimension is still lacking.
In this article, a new temperature-dependent heat source is developed to simulate deep penetration in SLM. The parameters of the heat source are no longer fixed during the simulation. At the beginning of each time step, the parameter of laser penetration is adjusted according to the pre-defined criteria. Experience shows that the temperaturedependent heat source needs to be calibrated only once for each material. The FE model with the new heat source is employed to simulate the SLM AlSi10Mg and SLM Inconel 718.
The organization of paper is as follows: (i) Section 2 first expounds the FE model using the new proposed heat source; (ii) as illustrative examples, Section 3 shows the comparisons of simulation results and experimental measurements [32,33]; and (iii) conclusions and perspectives are presented in Section 4. Figure 2 shows the dimensions of the model and schematic diagram process of thermal behavior in the SLM process. Thanks to the symmetrical plane, only half of the model is used. The powder layer (orange) is considered as a continuum medium with effective thermal conductivity and density, and the powder state can transform to a compact state once the temperature is superior to temperature fusion. This powder-compact transition is achieved by phase transformation [34,35]. The heat losses due to convection and radiation are also be taken into account.

Governing Equations
The governing equation of heat transfer in FEM [36], coupled with the change of state (powder, compact material), can be expressed as where ρ = ∑ phases p i ρ i , K = ∑ phases p i K i and C = ∑ phases p i C i . p i ρ i K i C i are, respectively, the proportion of each state (powder, compact material); the density, conductivity, and specific heat of state i (i = 1 is powder state, i = 2 is compact state); T is the temperature at the nodes; and q is the function to describe the heat generated by the laser. Finally, the convective and radiative losses are taken into account by the following boundary conditions: where h conv , σ, are the coefficient exchange with air, Stefan-Boltzmann constant, and the surface emissivity, respectively. T and T air are in degrees Kelvin.

FE Models and Meshes
Figures 3 and 4 show the FE meshes for SLM AlSi10Mg and SLM Inconel 718, respectively. Mesh refinement only appears along the trajectory of the heat source, and the minimum element size in the substrate is 5(length) × 5(width) × 9(depth) µm.    [32] and Liu et al. [33] are used as references, therefore, the thicknesses of the powder layer are 20 µm for AlSi10Mg and 40 µm for Inconel 718. The apparent heat capacity method is employed for Inconel 718 to incorporate the latent heat (L = 210 kj * kg −1 ) due to the phase changing between solid and liquid [37]. The substrate material properties of Al2024 used in SLM AlSi10Mg are exactly the same as Liu et al. [33]. For the sake of convergence, the thermal conductivity of powder varies linearly from 20 degrees to temperature fusion. The thermal properties of powder at ambient temperature (20 degrees) are calculated from the technique presented by Loh et al. [11], and the effective properties of powder are described by the following equations: where Φ is the porosity of powder. In this study, values of Φ = 0.25 and n = 4 are taken [11]. Since the pure heat transfer model has neglected fluid flow in melt pool and multiphysics problems related to the keyhole (see Figure 6), which leads to an under-estimation of molten pool depth and over-estimation of molten pool width, in the literature, the anisotropic thermal conductivity technique has been reported to be useful to improve the precision of molten pool dimensions. This technique has been used to simulate deep penetration by multiplying thermal conductivity by an enhancement factor [33,38,39], as shown in Equation (5).
where K xx , K yy , K zz represent the thermal conductivity in each direction. λ xx , and λ yy and λ zz are the temperature-dependant enhancement factors of k (presented in Table 5).
The depending temperature (T c ) and the numerical value of (M) should be estimated and validated from experiments. The enhancement factors are defined according to the following equation: where T c is greater than temperature fusion and M is generally set to 2-30 [40]. The melt pool depth is supposed to be in the z direction. Liu et al. [33] reported that values of T c = 750 degrees and M = 15 give a good correlation with experiments in SLM AlSi10Mg; therefore, the same values are used in the present study. However, the enhanced factors are calculated by comparing the predicted values with the experimental results, which is quite complex and time consuming. Moreover, a sharp evolution of thermal conductivity with temperature in one direction would probably lead to a difficulty in convergence.
The principle of simulating deep penetration was inspired by the 3D volumic heat source. The model we proposed is capable of increasing the parameter of penetration automatically in the heat source equation to consider the neglected transport phenomena in the melt pool, and the energy re-distribution in the depth direction related to the keyhole is represented by introducing a parameter that defines the position of maximal power in the depth direction. More details are presented in Section 2.3.

Heat Source Modeling
Considering the interactions between the powder and laser, the laser penetration depth η is generally equal to the thickness of the powder layer. Therefore, the 3D Gaussian heat source is represented as where r 0 , λ, P, r, and z are, respectively, the laser beam radius, energy absorptivity, laser power, the distance between the beam center and the point, and the distance from the top surface.
In case of deep penetration, since the governing Equation (1) neglects the fluid flow within the melt pool and the effects of the keyhole (see Figure 6), the melt pool size simulated with the Gaussian heat source described by Equation (7) will therefore underestimate the penetration and over-estimate the width. We can start with a comprehensive analysis related to these neglected effects.

1.
The fluid flow within the melting pool (Figure 6a) can lead to a more homogeneous temperature distribution by transporting the received energy with the mass flow. Therefore, the anisotropic thermal conductivity technique is well suited to consider the transport effect. In case of high energy input, higher-temperature fluid from the surface heated by a laser is carried toward the bottom of the melt pool [41], which results in a deep-penetration region in the central melt pool. Thus, we propose to increase the parameter η to simulate the penetration. Since η evolves at each time step, η t is now used in Equation (8). 2. The effect of multiple reflections on the keyhole (Figure 6b) has been studied by Han et al. [42]. This effect changes the energy distribution in depth, and the maximal power is no longer present at the top surface, which can lead to deeper penetration. One should note that this effect cannot be represented by imposing an enhancement factor for conductivity. The energy distribution described by Equation (7) should be modified by replacing exp − |z| where z 0 represents the position of maximal power in the depth direction. This parameter is not a fixed value due to the instability of the keyhole, and it is impossible to measure this value from experiments at the moment. Therefore, we propose z 0 , which is proportional to η t , for simplicity.
To conclude, the thermal effects related to Figure 6a can be compensated by increasing the parameter of η t , while η t cannot be superior to the melt pool depth. The η t is calculated in an explicit manner, where η t is updated based on the average temperature of controlled volume (see Figure 7) calculated from instant t − dt. The effect of multiple reflection (Figure 6b) can be represented by employing exp − |z−z 0 | η t in Equation (8), and the melt pool depth will be equally enhanced.
In order to ensure a constant energy input as well as consider the variation of η t and z 0 , Equation (7) is transformed into the following form: q(x, y, z) = ; if z not greater than η 0 ; otherwise (8) where N is recalculated at the beginning of each time step to ensure a constant energy input, and η t is also recomputed at beginning of each time step according to the average temperature calculated from selected nodes. We now focus on the procedure of calculating η t , and this computation is achieved in two steps.

2.3.1.
Step 1: Node Selection Procedure r 0 represents the radius of the laser beam. r 1 is is set to 1 3 r 0 , since more than 50% of the energy is distributed within this small area. If we take the moment t 0 as example, the center and height of the heat source are (0, y t 0 , 0) and η t 0 . A control volume of cylindrical shape is defined by r 1 , h and the center of the cylinder (0, y t 0 , η t 0 ), where r 1 is the radius and h represents the height of cylinder (see Figure 7). We calculate an average temperature from all the nodes included in control volume, and this average temperature is used as the criterion for adjusting η t .  Figure 8 shows the general procedure for selecting the nodes. A program (SIL) will loop on all the nodes present on the laser path, and the selected nodes should satisfy two conditions: where (x, y, z) are the coordinates of the node. One should note that h/2 should be larger than the minimal element size in the depth direction; otherwise, no nodes could be selected for Step 2. However, h/2 should not be so large as to induce the instability due to the adjustment of η t . In this study, h is set to 20 µm. Finally, the green nodes shown in Figure 8 are automatically selected for step 2. Compared to the strategy of selecting only one node at the bottom center of the laser beam, selecting numerous nodes can provide more representative results, and the average temperature is supposed to be less sensitive to oscillations related to time steps.

Step 2: Adjusting the Heat Source Parameters
The average temperature (T average ) of selected nodes is calculated from the results of the previous time step. Since the parameter η t cannot be greater than the melting pool depth, the material melting temperature is used as a criterion. Figure 9 presents the possible case according to the relationship between T average and T melting . The possible cases at instant t are as follows: (a) If T average T melting , η t should be increased for next time step. If T average 2T melting , T melting T average T melting , η t+1 keeps the same value (η t ) for the next time step.
(c) For the case T average < 2 3 T melting , which means that η t+1 is over-estimated, η t+1 = Once a new η t+1 is determined, N should be calibrated again to ensure a constant energy input during the process.  Figure 10 shows the pseudocode for controlling the simulation, which is also a summary of the temperature-dependent heat source. The simulations are performed in two stages: (i) the heat source's parameters are adjusted in the transient stage, and a small time step is used; (ii) the heat source becomes stable in the steady-state stage, and a larger time step can be employed. Figure 11 provides an overview of the flowchart of the proposed modeling flow of the SLM process. One should note that steady-state conditions require the same heat source parameters in N successive transient stages to be obtained, and we take N = 5 in this study. In this article, SYSWELD TM software and its Marco language SIL are used to implement the proposed heat source. The same procedure can also be achieved using other software and Macro languages, such as Python for Abaqus and APDL for Ansys.

Results and Discussions
The melt pool width and depth are illustrated in Figure 12, and both experimental results and numerical results are measured in the same way. To ensure the good stability of prediction during the transient stage, a time step of smaller than r 0 3 * V scanning is recommended. FE software SYSWELD [43] was utilized to model SLM.  Table 1 presents the scanning speed and laser power applied in SLM simulations. The other parameters, such as emissivity, coefficient exchange, initial conditions, and so on are the same as those used by Liu et al. [33]. The experimental results of sample 1# are used to calibrate z 0 (z 0 = 0.66 * η t ). The simulation results of V200-P200 (scanning speed labeled V: 200 mm/s and laser power labeled P: 200 Watt) are used to illustrate certain characters of the proposed heat source. Initial temperature ( • C) 25 Figure 13 shows the melting pool shape at different moments during the simulation, and the contours are shown at the symmetrical plane. The melting pool depth increases progressively over time, since the heat source increases its penetration automatically.

SLM AlSi10Mg
The final simulated melt pool shape is shown in Figure 14, and two stages can be clearly identified. The melt pool depth keeps increasing in the transient stage, while the melt pool depth is almost constant in the steady-state stage. A slight increase of the melt pool depth can be observed in the steady-state stage, which is perhaps due to the boundary effect as the model is short.  Another investigation concerns the sensibility of the time step, as it is known that η t is calculated in an explicit way. Figure 15 shows the molten profiles simulated with different time steps, where the maximal recommended time step is dt = r 0 3 * V scanning = 0.06 ms.
The comparison of the melt pool depths shows that a temperature-dependent heat source is stable and slightly sensitive to the time step. A maximum difference of 3.8% is observed, and the increase of melt pool depth in dt = 0.125 ms is perhaps due to the boundary effect. Even though the time step is different in this study, the number of time steps required to achieve the steady-state stage is almost the same (about 29 transient time steps).  Figure 16 provides the comparison of the cross-section morphology of three melt pool shapes. The dimensions of the melt pool given by experiment and simulation are quite similar, and the relative error between the experimental and simulated data is 4.8% for width and 2% for depth. One should note that the heat source parameters are adjusted automatically without any manual intervention. The comparison of the molten pool dimensions of the experimental and numerical results is presented in Figure 17, and good correlations among experimental results (Exp-Width and Exp-Depth) [33], numerical results with the proposed heat source (Auto-Width and Auto-Depth), and with the anisotropic conductivity method (Anisotropic-Width and Anisotropic-Depth) [33] are found. The proposed heat source is capable of predicting the melting pool dimension with different process parameters correctly. Compared to the experimental results, the average relative errors given by the proposed model are generally less than 4.8% for depth and 5.2% for width. A better precision of prediction than the anisotropic conductivity method is thus observed. Moreover, with fixed laser power, both the width and depth of the molten pool decrease if the scanning speed increases. When the scanning speed is fixed at 600 mm/s, the width and depth increase with increasing laser power.

SLM Inconel 718
We now apply the proposed heat source for the SLM Inconel 718, and the experimental results are taken from Sadowski et al. [32]. The initial parameter η t=0 is equal to 40 µm, and z 0 is 0.4 * η t . The reported absorptivity of Inconel 718 varies from 0.3-0.87 [27,44]; therefore, 0.7 is used in this study. The processing parameters (scanning speed and laser power) used are shown in Table 2, and other parameters can be found in the work of Romano et al. [27]. Initial temperature ( • C) 20 Figure 18 presents the melting pool at different moments during the simulation. Compared to the temperature distribution in Figure 13, the heat diffuses very slowly in Inconel 718 because of its low diffusivity. Figure 19 shows the final simulated molten pool shape (speed: 200 mm/s, power: 100 Watt). The molten pool depth increases in the transient stage and keeps constant in the steady-state stage. Boundary effects become negligible due to low diffusivity. Thus, the melt pool depth in the steady-state stage does not increase as we observed in SLM ALSi10Mg.
Similarly, the effect of the time step is studied in SLM Inconel 718. Figure 20 presents the molten profiles after the SLM process and the variation of melt pool depth in the scanning direction. In the steady-state stage, a maximum difference of 2.9% is observed, and we can conclude that the precision is satisfactory.
Finally, a comparison of the molten pool dimensions given by experiments and numerical simulations is presented in Figure 21. The simulated melt pool depth gives good precision compared to experiments, while the width is slightly under-estimated. Thus, a maximal relative error of 9.6% is found.    The proposed heat source is therefore validated through the comparisons of its results (melt pool morphology for two materials: AlSi10Mg and Inconel 718) with those of the existing numerical method (anisotropic thermal conductivity technique) and experimental tests. Compared to the anisotropic thermal conductivity technique, the procedure of parameter calibration in the proposed heat source (only the parameter z 0 must be calibrated) is easier. Moreover, numerical investigation shows that the proposed heat source is slightly sensitive to the choice of the time step. In contrast, a sharp increase of thermal conductivity induced by the anisotropic thermal conductivity technique leads to the difficulty of converging, which results in only very small time steps being applied. The difficulty of converging could be further increased if the metallurgical transformation is considered. Therefore, the proposed heat source is more efficient in terms of computation cost than the anisotropic thermal conductivity technique. Note that the proposed heat source can be subsequently used for other materials and processes (for example, laser welding by the keyhole technique).

Conclusions
A thermal model based on a new heat source has been developed to simulate deep penetration during the SLM process. This heat source can adjust its parameters during the simulation to compensate the neglected thermal effects related to the fluid flow and keyhole. Its parameters become fixed once the temperature distribution becomes stable. The proposed model has been validated on experimental results (for two materials: ALSI10Mg and Inconel 718) from the literature. The comparisons of melt pool morphology showed a good correlation between numerical and experimental results.

Conflicts of Interest:
The authors declare no conflict of interest.