Transient Powder Melting in SLM Using an Analytical Model with Phase Change and Spherical Symmetry in a Semi-Infinite Medium

In this work, we introduce an analytical expression for approximating the transient melting radius during powder melting in Selective Laser Melting (SLM) assumed with a stationary laser heat source. The purpose of this work is to evaluate the suggested analytical approach in determining the melt pool geometry during laser processing, by considering heat transfer and phase change effects. This will allow for the rendering of the first findings on the way to a quasi-real time calculation of the melt pool during laser melting, which will contribute significantly to the process design and control, especially when new powders are applied. Initially, we consider the heat transfer process associated with a point heat source, releasing a continuous and constant power (in a semi-infinite powder bed. On the point of the heat source the temperature is infinite, and the material starts to melt spherically outwards, creating an interface that separates the solid from the molten material; we assume different properties between the two phases. Unlike the cases of the cartesian and cylindrical coordinates, (in a cartesian coordinate the heat source is over a plane, i.e., W/m2, and in cylindrical along a line, i.e., W/m), where the melting process is proportional to the square root of time, in spherical coordinates the melting stops at a finite radius, i.e., a maximum radius, which depends only on the heat source, the conductivity of the solid and the difference between the far-field temperature and the melting temperature of the material. Here we should also point out that to achieve continuous melting in spherical coordinates the power of the source must increase with the square root of the time. The obtained analytical expression for the maximum melting radius and the approximate expression for its dependence on the time compare well with the numerical results obtained by a finite element analysis.


Introduction
Melting due to a heat source is very significant in thermal manufacturing processes such as Selective Laser Melting (SLM) [1][2][3][4][5][6] and Selective Laser Sintering (SLS) [7,8], spot welding, torch welding [9][10][11], and arc welding [12], to name a few. In these processes, an understanding of the melting process is very important due to the many other physical phenomena taking place simultaneously, such as Marangoni effects and the internal flow in the molten material.
The fundamentals of the thermodynamic change phase of solids are presented comprehensively by Tosun [13], whereas in the work of Berveiller & Fischer [14] a deeper insight into the influence of the change phase on the mechanical properties of solids is outlined. Here, kinematic models to describe the phase change effects in metals during melting and solidification, especially in the case 2 of 9 of steel alloys, are expedient [15]. During the phase change from solid to liquid and back to solid, thermal and plastic stresses occur due to a lower yield strength in metals at higher temperatures [16]. After a cooling down to an ambient temperature due to a non-linear interaction of thermal, elastic and plastic stresses and the phase change, rest stresses are formed in the solid structure, so called residual stresses [17,18]. In particular in the case of SLM processes, a residual stresses formation is of great importance since residual stresses are dominated by the actual process parameters [19]. Modeling methods are applied in order to simulate the formation of residual stresses that have a significant impact in process feasibility and product properties [20,21]. A deeper insight into the characterization and evaluation of residual stresses by means of experimental measurements and analytical models is presented by Ghidelli et al. [22].
Fundamental work on phase change dates back to Lamè & Clapeyron [23] and to Stefan [24] who consider the problem of ice formation. Phase change in cylindrical and spherical coordinates were later considered by Frank [25], Paterson [26], and Cho & Sunderland [27]. Paterson [26] addressed the problem of a line heat source in cylindrical coordinates and obtained an expression for the propagation of the melting interface. The temperature along the line heat source is infinite, and an interface separates the material into two regions: the molten phase and the solid phase. The velocity of the interface is proportional to the square root of time, and the proportionality constant is obtained through a characteristic equation by employing the boundary condition describing the energy conservation along the boundary of fusion. The analytical expression for the process was obtained by employing a similarity variable that reduced the partial differential equations (PDE) into ordinary differential equations (ODE). Finally, the characteristic algebraic equation is obtained.
Analytical expressions for the phase change process can be also obtained in one-dimensional cartesian coordinates in the case of the solidification of a supercooled liquid, and the melting or solidification in a half-space with a constant temperature along the boundary [28]. As mentioned before, a common characteristic of many one-dimensional problems, particularly in cartesian coordinates, is the reduction of the PDE to an ODE using the characteristic variable x/ √ t. Detailed accounts of many problems with different configurations in solidification and melting are described in Alexiades & Solomon [29], Carslaw & Jaeger [30], Ghez [31] and Hu & Argyropoulous [32], along with the methodologies for addressing more complicated problems, both analytically and approximately.
An analytical solution for solidification/melting in spherical coordinates is not available. The physical reason for this is that the temperature drops rapidly enough along the radius, and after a certain point the solidification/melting cannot be sustained due to the reduction of the temperature gradient. Mathematically, the situation can be easily understood by considering the problem of one dimensional heat conduction in an infinite medium in the three basic coordinate systems, i.e., cartesian, cylindrical and spherical. Although for steady-state conditions there is an analytical solution in spherical coordinates and the temperature is proportional to the inverse radius 1/r, this is not the case in cylindrical and cartesian coordinates; in cylindrical coordinates, the temperature is proportional to the log[r] and in artesian coordinates it is proportional to x [33], i.e., for cartesian and cylindrical coordinates the only possible solution is time dependent. Hence Paterson [26], who considered the problems of phase change in cylindrical and spherical coordinates, and was able to obtain an analytical solution in spherical coordinates only under the assumption that the power of the heat source increases with the square root of time; hence, this solution is of less practical importance.
In the existing literature, time-intensive thermal models of the transient temperature distribution during SLM processing are presented using equivalent heat sources on the basis of the finite element analysis [34,35]. In this work, we focus exclusively on the transient heat transfer effects during the phase change, i.e., the melting, of solid powder to liquid, by proposing an analytical model. We consider the problem of phase change in spherical coordinates assuming a continuous but constant heat source at a single point. This simplifies the problem, as the solution is spherically symmetric. In what follows, we present the problem, the steady-state solution and a transient approximate analytical solution, and compare them with a numerical simulation obtained using a finite element analysis. Here, the main goal is to appraise the proposed analytical solution, which provides a quasi-real time estimation of the melt pool. Future enhancements are intended to assist the SLM process, as well as the product design and control, especially due to the recent rising industrial application of additive manufacturing (AM) technologies [36].

Problem Statement and Solution
Consider the case of semi-infinite solid material, i.e., powder beds, that we wish to melt in order to subsequently produce a solid component. A point laser beam provides a heat source of power . Q that melts the material with a spherical symmetry. The configuration is shown in Figure 1. We assume that the solid and the molten material have the same density, so the shape of the configuration remains the same, i.e., no annulus is created. However, we assume that the heat capacity and the conductivity of the two phases are different. This leads to a discontinuity of the temperature gradient at the interface (R[t]) between the solid and the liquid (molten material). The interface (R[t]) separates the solid from the molten material and is a function of time. real time estimation of the melt pool. Future enhancements are intended to assist the SLM process, as well as the product design and control, especially due to the recent rising industrial application of additive manufacturing (AM) technologies [36].

Problem Statement and Solution
Consider the case of semi-infinite solid material, i.e., powder beds, that we wish to melt in order to subsequently produce a solid component. A point laser beam provides a heat source of power Q  that melts the material with a spherical symmetry. The configuration is shown in Figure 1. We assume that the solid and the molten material have the same density, so the shape of the configuration remains the same, i.e., no annulus is created. However, we assume that the heat capacity and the conductivity of the two phases are different. This leads to a discontinuity of the temperature gradient at the interface (R[t]) between the solid and the liquid (molten material). The interface (R[t]) separates the solid from the molten material and is a function of time. The equations that describe the two phase-problem are as follows [28]: (i) The liquid phase: (ii) The solid phase: The equations that describe the two phase-problem are as follows [28]: (i) The liquid phase: 4 of 9 (ii) The solid phase: (iii) Energy balance on the interface: where α is the thermal diffusivity, k is the thermal conductivity, L is the latent heat, ρ is the density that is constant, i.e., ρ s = ρ l , and T melt is the melting temperature of the material. In the above equations, the subscript l denotes the liquid phase (molten material) and s the solid (powder).
The above system of equations and boundary conditions (Equations (1)- (7)) have no analytical solution. There is only an analytical/similarity solution if the heat source increases with the square root of time [26], i.e., However, if we assume steady-state conditions, i.e., we drop the time derivative terms from the system of Equations (1)-(7), and we obtain the following steady-state solutions: T l , T s , and R for the temperature fields and the radius, respectively: and It is important to point out that the above steady-state solutions must also be valid in the case of different densities between the solid and the liquid. This is justified because at a steady-state the terms related to convection must be zero, and the system of Equations (1)-(7) would be applicable.

Approximate Analytical Solution
An approximate solution for the interface R[t] can be obtained by considering that the process is controlled by diffusion in the solid phase. This is justified because (i) the powder has a smaller thermal diffusivity than the molten material, (ii) the melt pool is confined in a small radius and (iii) the steady-state solution of the interface R (Equation (10)) is independent of the properties of the liquid. Hence, for the case of a constant, continuous, thermal energy source . Q at the origin of a semi-infinite, homogeneous medium, the temperature distribution is obtained as [37]: Hence, the location of the interface R[t] with the temperature T = T melt is given by numerically solving the implicit equation:

Numerical Solution
The finite element simulation was performed on the commercial software ANSYS Workbench [38]. For the properties of the solid material, we have used the properties of the powder IN718 [2], i.e., k s = 0.37 W/(mK), α s = 2.76 × 10 −7 m/s 2 , T ∞ = 20 • C and T melt = 1300 • C. For a heat source of power . Q = 4 W, using Equation (10), we obtain that at a steady-state the interface is located at R = 0.13 mm.
In the FE simulation, the heat power Hence, the location of the interface R[t] with the temperature T = Tmelt is given by numerically solving the implicit equation:

Numerical Solution
The finite element simulation was performed on the commercial software ANSYS Workbench [38]. For the properties of the solid material, we have used the properties of the powder IN718 [2], i.e., ks = 0.37 W/(mK), αs = 2.76 × 10 −7 m/s 2 ,  T = 20 °C and Tmelt = 1300 °C. For a heat source of power Q  = 4 W, using Equation (10), we obtain that at a steady-state the interface is located at R = 0.13 mm.
In the FE simulation, the heat power Q  was defined on an infinitesimal surface of radius 0.

Numerical Results with no Phase Change
To verify the numerical analysis, we first addressed the problem of the continuous release of heat at single point on a homogeneous semi-infinite medium, i.e., without considering the varying material properties at the melt pool interface (Equation (7)). We have used the properties of the solid mentioned earlier. The analytical solution for this problem is given by Equation (11), hence we expect the numerical solution to match Equation (12). In Figure 3, we compare the results of the simulation with Equation (12). The figure shows the location of the radius where the temperature is at T melt for the two cases, i.e., numerical (dashed curve) and analytical (solid curve, Equation (12)). For the given data, the steady-state radius is R = 0.00134 m. As expected, the results are very accurate, and the error is less than 1%.
To verify the numerical analysis, we first addressed the problem of the continuous release of heat at single point on a homogeneous semi-infinite medium, i.e., without considering the varying material properties at the melt pool interface (Equation (7)). We have used the properties of the solid mentioned earlier. The analytical solution for this problem is given by Equation (11), hence we expect the numerical solution to match Equation (12). In Figure 3, we compare the results of the simulation with Equation (12). The figure shows the location of the radius where the temperature is at Tmelt for the two cases, i.e., numerical (dashed curve) and analytical (solid curve, Equation (12)). For the given data, the steady-state radius is R = 0.00134 m. As expected, the results are very accurate, and the error is less than 1%.  (12)), and the dashed curve is the numerical solution obtained using the finite element analysis (ANSYS [38]).

Numerical Results with Phase Change
The numerical analysis is enhanced with the definition of the varying, i.e., non-linear specific heat capacity between the solid and liquid phases, implying that the phase change from the solid to liquid state requires the latent heat at the melting temperature. We have assumed a pure material, i.e., the phase change occurs at the distinct temperature Tmelt. The liquid phase and the solid phase are characterized by different properties of the heat capacity, specifically cps = 351 J/(kg·K) and cpl = 643 J/(kg·K), resulting in a latent heat of 4.59 × 10 5 J/kg. In addition, the density of the liquid is ρl = 7756 kg/m 3 and kl = 26.63 W/(m·K). A parameter sensitivity analysis shows that the thermal conductivity variation due to the phase change does not influence the melt pool creation significantly, i.e., the dominant material property in this case is the heat capacity, as justified by literature findings [39]. In Figure 4, we show the results of the numerical simulation, and we compare them with the approximate transient analytical result (Equation (12)). It is evident that both the numerical and approximate results (Equation (12)) approach the steady-state solution (Equation (10)), the only difference being the rate at which this is approached. It is interesting, however, that although the approximate solution (Equation (12)) assumes a single material, in particular solid powder, it provides results close to the numerical solution, the maximum error being approximately 10%. This is not surprising as the melt pool is isolated in a small region.
In addition, because of the latent heat and the higher heat capacity, the interface obtained by the numerical solution requires more time to achieve a certain radius, hence the radius obtained from the analytical solution is always higher than the one obtained from the numerical solution.

Numerical Results with Phase Change
The numerical analysis is enhanced with the definition of the varying, i.e., non-linear specific heat capacity between the solid and liquid phases, implying that the phase change from the solid to liquid state requires the latent heat at the melting temperature. We have assumed a pure material, i.e., the phase change occurs at the distinct temperature T melt . The liquid phase and the solid phase are characterized by different properties of the heat capacity, specifically cp s = 351 J/(kg·K) and cp l = 643 J/(kg·K), resulting in a latent heat of 4.59 × 10 5 J/kg. In addition, the density of the liquid is ρ l = 7756 kg/m 3 and k l = 26.63 W/(m·K). A parameter sensitivity analysis shows that the thermal conductivity variation due to the phase change does not influence the melt pool creation significantly, i.e., the dominant material property in this case is the heat capacity, as justified by literature findings [39]. In Figure 4, we show the results of the numerical simulation, and we compare them with the approximate transient analytical result (Equation (12)). It is evident that both the numerical and approximate results (Equation (12)) approach the steady-state solution (Equation (10)), the only difference being the rate at which this is approached. It is interesting, however, that although the approximate solution (Equation (12)) assumes a single material, in particular solid powder, it provides results close to the numerical solution, the maximum error being approximately 10%. This is not surprising as the melt pool is isolated in a small region.
In addition, because of the latent heat and the higher heat capacity, the interface obtained by the numerical solution requires more time to achieve a certain radius, hence the radius obtained from the analytical solution is always higher than the one obtained from the numerical solution.  (12)), and the dashed curve is the numerical solution. In the latter, we assume that there are two phases: solid and liquid.

Conclusions
The findings of our work are summarized as follows: (1) We consider the heat conduction problem associated with a continuous source of power at a single point in a semi-infinite material. At the point of the heat source, the temperature is infinite, and the melting process is spherically symmetric. Unlike cartesian and cylindrical coordinates, there is no analytical solution; there is, however, a steady-state solution, i.e., the melting process reaches a maximum radius. The radius is proportional to the power of the heat source and inversely proportional to the conductivity of the solid and the difference between the melting temperature and the temperature at infinity, as proposed in the work of [26]. (2) An approximate analytical solution of the melting radius as a function of time is obtained by assuming a single material, i.e., solid powder, and by locating the radius where the temperature is at the melting temperature, as in [37]. (3) This simple approximation has the same steady-state result as the numerical solution, and also provides transient results close to the numerical solution, although the numerical solution includes the latent heat, and a higher heat capacity and conductivity for the fluid similar to the finite element analyses conducted by [34,35] for a moving heat source (Figure 4). The reason for this is that the heat conduction process is controlled by the material with the lower thermal diffusivity, i.e., the solid powder, and that the melt pool has small dimensions. The difference between the two is that the numerical result requires more time to achieve the steady-state solution because of the extra energy required due to the latent heat and the non-linear heat capacity of the molten material.  where the temperature has the constant value T melt . The solid curve is the analytical solution obtained assuming a single material, i.e., homogeneous powder (Equation (12)), and the dashed curve is the numerical solution.
In the latter, we assume that there are two phases: solid and liquid.

Conclusions
The findings of our work are summarized as follows: (1) We consider the heat conduction problem associated with a continuous source of power at a single point in a semi-infinite material. At the point of the heat source, the temperature is infinite, and the melting process is spherically symmetric. Unlike cartesian and cylindrical coordinates, there is no analytical solution; there is, however, a steady-state solution, i.e., the melting process reaches a maximum radius. The radius is proportional to the power of the heat source and inversely proportional to the conductivity of the solid and the difference between the melting temperature and the temperature at infinity, as proposed in the work of [26]. (2) An approximate analytical solution of the melting radius as a function of time is obtained by assuming a single material, i.e., solid powder, and by locating the radius where the temperature is at the melting temperature, as in [37]. (3) This simple approximation has the same steady-state result as the numerical solution, and also provides transient results close to the numerical solution, although the numerical solution includes the latent heat, and a higher heat capacity and conductivity for the fluid similar to the finite element analyses conducted by [34,35] for a moving heat source (Figure 4). The reason for this is that the heat conduction process is controlled by the material with the lower thermal diffusivity, i.e., the solid powder, and that the melt pool has small dimensions. The difference between the two is that the numerical result requires more time to achieve the steady-state solution because of the extra energy required due to the latent heat and the non-linear heat capacity of the molten material.