Analytical Thermal Modeling of Metal Additive Manufacturing by Heat Sink Solution

Metal additive manufacturing can produce geometrically complex parts with effective cost. The high thermal gradients due to the repeatedly rapid heat and solidification cause defects in the produced parts, such as cracks, porosity, undesired residual stress, and part distortion. Different techniques were employed for temperature investigation. Experimental measurement and finite element method-based numerical models are limited by the restricted accessibility and expensive computational cost, respectively. The available physics-based analytical model has promising short computational efficiency without resorting to finite element method or any iteration-based simulations. However, the heat transfer boundary condition cannot be considered without the involvement of finite element method or iteration-based simulations, which significantly reduces the computational efficiency, and thus the usefulness of the developed model. This work presents an explicit and closed-form solution, namely heat sink solution, to consider the heat transfer boundary condition. The heat sink solution was developed from the moving point heat source solution based on heat transfer of convection and radiation. The part boundary is mathematically discretized into many heats sinks due to the non-uniform temperature distribution, which causes non-uniform heat loss. The temperature profiles, thermal gradients, and temperature-affected material properties are calculated and presented. Good agreements were observed upon validation against experimental molten pool measurements.


Introduction
Powder bed metal additive manufacturing (PBMAM) can produce geometrically complex parts with effective cost. For instance, with the use of powder bed metal additive manufacturing (PBMAM) configuration, high-density laser power is employed to fully melt and fuse metal powders to build parts in a layer-by-layer manner. The high thermal gradient due to the repeated rapid heating and solidification cause defects in the produced parts, such as cracks [1], porosity [2,3], undesired residual stress [4,5], and part distortion [6]. Different techniques have been developed to monitor and control temperature conditions, namely experimental measurement, finite element method (FEM)-based numerical modeling, and physics-based analytical modeling.
In situ temperature measurements provide real-time temperature measurements during the heating state and cooling state in different additive manufacturing (AM) processes. Embedded thermocouples have commonly been used to measure the temperature on or inside the substrate [7]. Thermal imaging cameras (infrared cameras) were used in previous work to measure the temperature This work presents an explicit and closed-form solution, namely heat sink solution, for heat transfer boundary conditions in analytical thermal modeling of PBMAM. The heat sink solution was developed based on the heat transfer mechanism for convection and radiation. The temperature distribution was predicted using the heat sink solution and the moving point heat source solution in PBMAM of Ti6Al4V. Different numbers of heat sinks were used in the temperature prediction to investigate the influence of heat sinks. The thermal gradient and material property variations were investigated based on the temperature prediction. The optimal number of heat sinks was determined by comparison with experimental measurements. With the same number of heat sinks, the temperature distributions were predicted and validated under various process conditions.

Materials and Methods
In this work, a closed-form solution, namely the heat sink solution, is presented to characterize the heat loss from the part boundary due to convection and radiation in PBMAM, as illustrated in Figure 1. The red arrow and green arrows represent heat input form the laser source and heat loss from convection and radiation at the part top boundary, respectively. The heat loss from the side boundary (x-z planes and y-z planes) is assumed to be negligible because of the significantly lower temperatures than the top boundary (x-y plane). L, W, and D denote the molten pool length, molten pool width, and molten depth at the laser heat source location, respectively. developed based on the heat transfer mechanism for convection and radiation. The temperature distribution was predicted using the heat sink solution and the moving point heat source solution in PBMAM of Ti6Al4V. Different numbers of heat sinks were used in the temperature prediction to investigate the influence of heat sinks. The thermal gradient and material property variations were investigated based on the temperature prediction. The optimal number of heat sinks was determined by comparison with experimental measurements. With the same number of heat sinks, the temperature distributions were predicted and validated under various process conditions.

Materials and Methods
In this work, a closed-form solution, namely the heat sink solution, is presented to characterize the heat loss from the part boundary due to convection and radiation in PBMAM, as illustrated in Figure 1. The red arrow and green arrows represent heat input form the laser source and heat loss from convection and radiation at the part top boundary, respectively. The heat loss from the side boundary (x-z planes and y-z planes) is assumed to be negligible because of the significantly lower temperatures than the top boundary (x-y plane). L, W, and D denote the molten pool length, molten pool width, and molten depth at the laser heat source location, respectively. Schematic drawing of the heat transfer mechanism in PBMAM. The red arrow represents heat input from the laser source. The green arrows represent the heat loss from part boundary due to convection and radiation. Here, x, y, and z denote coordinate directions, while L, W, and D denote the molten pool length, width, and depth, respectively.
The heat balance equation was used as the governing equation for the derivations of the heat source solution and heat sink solution. This equation describes the change of temperatures due to the energy input from a moving volumetric heat source, and thus can be employed for the PBMAM [31]. It is expressed as the following where u is internal energy, H is enthalpy, ρ is density, k is conductivity, is volumetric heat source, t is time, x is the distance from the heat source, V is the heat source moving velocity along the xdirection, and T is temperature. The heat conduction equation was derived from the heat balance equation with = 0, and = c as = • ) + , which can be rewritten by introducing the thermal diffusivity ( = /ρc) as the following Figure 1. Schematic drawing of the heat transfer mechanism in PBMAM. The red arrow represents heat input from the laser source. The green arrows represent the heat loss from part boundary due to convection and radiation. Here, x, y, and z denote coordinate directions, while L, W, and D denote the molten pool length, width, and depth, respectively.
The heat balance equation was used as the governing equation for the derivations of the heat source solution and heat sink solution. This equation describes the change of temperatures due to the energy input from a moving volumetric heat source, and thus can be employed for the PBMAM [31]. It is expressed as the following ∂ρu ∂t where u is internal energy, H is enthalpy, ρ is density, k is conductivity, . q is volumetric heat source, t is time, x is the distance from the heat source, V is the heat source moving velocity along the x-direction, and T is temperature. q, which can be rewritten by introducing the thermal diffusivity (κ = k/ρc) as the following where κ is thermal diffusivity (κ = k/ρc, where k, ρ, c are thermal conductivity, density, and specific heat, respectively), and x, y, z are the distances from the heat source. The transient-state moving point heat source solution was derived from the heat conduction equation by Carslaw and Jaeger [31]. The laser heat source was assumed as a moving point heat source for a semi-infinite medium. It is expressed as the following where P is heat source power, η is absorption, t is the current time, t is previous time, and x, y, z are the corresponding distances from the laser source. The transient-state moving point heat source solution can be rewritten by integrating t' from 0 to t as the following where R 2 = x 2 + y 2 + z 2 is the total distance from the laser source, ξ is a time-related integration variable, which are introduced for concise expression and easy implementation of the heat source solution.
The steady-state moving point heat source solution can be derived with infinite t. It is expressed as The heat sink solution was derived with equivalent power for heat loss from convection and radiation, and zero moving velocity. The convection and radiation can be calculated as where Q conv and Q rad denote heat loss due to convection and radiation, respectively, A is the area of the heat sink, h is the convection coefficient, ε is emissivity, σ is the Stefan-Boltzmann constant, T is the temperature of the heat sink that can be estimated by the moving point heat source solution, and T 0 is room temperature. The equivalent power for heat loss at the part boundary can be expressed as Each heat sink is a portion of the part boundary that does not move as the laser heat source. Therefore, the moving velocity of each heat sink becomes zero (V = 0). The heat sink solution is expressed as Materials 2019, 12, 2568

of 15
The part boundary is mathematically discretized into many sections (heat sinks), considering the non-uniform temperature distribution at the part boundary, which causes non-uniform heat loss at the boundary.
where n is the number of heat sinks, which can be determined with experimental calibration to accurately calculate the heat loss at the part boundary, and i is the index of each heat sink.
The final temperature solution is constructed from the superposition of heat source solution and heat sink solutions as the following In addition, the latent heat (L f ) is considered using the heat integration method [30], in which the temperatures are lowered because phase transformation takes place with continuous heat input.
The presented heat sink solution provides a computationally efficient method to consider the heat transfer boundary conditions in PBMAM without resorting to FEM or any iteration-based simulations. The analytical model is constructed from the superposition of closed-form solutions, namely the heat sink solution and moving point heat source solution. The promising short computational time allows temperature prediction for larger scale parts and process-parameter planning through inverse analysis. Therefore, the developed model has improved usefulness in real applications.

Results and Discussion
In this work, the three-dimensional temperature distribution was predicted by the presented model in powder bed metal additive manufacturing (PBMAM) of Ti6Al4V. The heat sink solution was derived and employed to impose the heat transfer boundary condition without significantly compensating the computational efficiency in analytical temperature modeling. The implementation of the heat sink solution was investigated by applying different numbers of heat sinks in the temperature prediction. The temperature profile, thermal gradient, and the temperature-affected material property variation were calculated in the single-track scan. The molten pool dimensions were calculated by comparing the predicted temperature to the material melting temperature. They were used to determine the optimal number of sinks in comparison with the documented experimental values. The computational time was recorded and presented. The molten pool dimensions were then calculated by the presented model with the optimal number of heat sinks under different process conditions. Validation of the calculated molten pool dimensions was included.
The presented model was constructed from the superposition of the heat source solution and heat sink solution. The temperature rise due to the heat input from the moving laser heat source was calculated using the moving point heat source solution. The power absorptivity in the heat source solution was adopted as 0.77 [12], which is related to the laser wavelength, laser-workpiece offset distance, powder material properties, and powder-packing-related surface roughness. The temperature drop due to the heat loss from the part boundary was calculated using the heat sink solution. The part top boundary was mathematically discretized into many sinks considering the non-uniform temperature distribution at the boundary, which led to non-uniform heat loss at the boundary according to Equations (6) and (7). The number of heat sinks was determined by calibration based on the experimental measurement of molten pool dimensions in the literature [2]. The documented molten pool width and depth were measured based on the solidification microstructure using an optical  Table 1. The process parameters in PBMAM of Ti6Al4V are given in Table 2. Different laser powers and scanning velocities were applied to the PBMAM. The details of the process parameters and documented molten pool dimensions are given in Table A1 in Appendix A. Table 1. Materials properties and thermophysical properties of Ti6Al4V [12,[42][43][44].

Name Symbol Value Unit
Density ρ 4428 kg/m 3 Thermal conductivity (powder at T 0 ) k p 6.6 W/(m·K) Note: Solid thermal conductivity and solid specific heat are only used in the investigation of materials' property variations. The heat sink solution was derived from the moving point heat source solution, which has a singularity issue in the heat source solution. In other words, the temperature becomes infinite in the heat source solution with zero distance (Equation (5)). Therefore, the heat sink solution has an inherent singularity issue. In order to properly implement the heat sink solution, the influence of heat sinks on the temperature prediction was fully investigated. The temperature profiles and temperature gradient were calculated in single-track scans under test 1 condition (P = 100 W; V = 500 mm/s). The numbers of heat sinks on the top boundary were chosen as 4 × 4, 6 × 6, 8 × 8, and 10 × 10, respectively. The area of the heat sink under each setting is identical to the total area/number of heat sinks. For example, the area of each heat sink with 4 × 4 setting is identical to the total area/16. The area of each heat sink with 8 × 8 setting is identical to the total area/64. The heat loss from the side boundary was not considered because of the significantly lower temperature compared to the top boundary near the heat source location (x = 0.8 mm, y = 0.5 mm). The three-dimensional temperature profiles are illustrated in Figure 2. As shown in the temperature profiles on the top boundary at the x-y plane, the predicted heat-affected zones considering the heat loss (solid lines) were smaller than those predicted without considering the heat loss (dashed lines) at each temperature level. The more heat sinks, the smaller the heat affected zone, and vice versa. The employed heat sink solution can reduce the overestimation of temperature levels, and thus improve the prediction accuracy. With the heat sink solution, the complete understanding of the heat transfer mechanism in PBMAM can be implemented conveniently and with computational efficiently. Therefore, the usefulness of the developed analytical model can be significantly improved in real applications. predicted heat-affected zones considering the heat loss (solid lines) were smaller than those predicted without considering the heat loss (dashed lines) at each temperature level. The more heat sinks, the smaller the heat affected zone, and vice versa. The employed heat sink solution can reduce the overestimation of temperature levels, and thus improve the prediction accuracy. With the heat sink solution, the complete understanding of the heat transfer mechanism in PBMAM can be implemented conveniently and with computational efficiently. Therefore, the usefulness of the developed analytical model can be significantly improved in real applications. The temperature gradient profiles were plotted in Figure 3 with different numbers of heat sinks under test 1 condition ( = 100 W; = 500 mm/s). The large temperature gradient was observed at the near heat source location (x = 0.8 mm y = 0.5 mm) and near heat sink location (marked as red circles). The materials property variations, namely the thermal conductivity and the specific heat, The temperature gradient profiles were plotted in Figure 3 with different numbers of heat sinks under test 1 condition (P = 100 W; V = 500 mm/s). The large temperature gradient was observed at the near heat source location (x = 0.8 mm y = 0.5 mm) and near heat sink location (marked as red circles). The materials property variations, namely the thermal conductivity and the specific heat, were plotted in Figures 4 and 5 respectively. The material property variation was caused by the temperature-dependent nature and the temperature variation. were plotted in Figures 4 and 5 respectively. The material property variation was caused by the temperature-dependent nature and the temperature variation.      To determine the proper number of heat sinks, the calculated molten pool dimensions were compared to the experimental values documented in the literature [2]. The molten pool dimensions were calculated by comparing the predicted temperature to the material melting temperature, as illustrated in Figure 6. The documented values were experimentally measured based on the solidification microstructure, as illustrated in Figure 7. To determine the proper number of heat sinks, the calculated molten pool dimensions were compared to the experimental values documented in the literature [2]. The molten pool dimensions were calculated by comparing the predicted temperature to the material melting temperature, as illustrated in Figure 6. The documented values were experimentally measured based on the solidification microstructure, as illustrated in Figure 7.     [2]. W and D denotes molten pool width and molten pool depth, respectively.
As shown in Figure 8, the closest agreement was observed with 6 × 6 heat sinks under test 1 conditions. The horizontal axis 0 × 0 denotes the calculated molten pool dimensions without a heat sink. In other words, the heat loss from the part boundary was not considered. Expt. denotes the experimental molten pool dimensions. The more heat sinks, the smaller the molten pool dimensions, Figure 7. Experimental measurements of molten pool dimensions based on the solidification microstructure [2]. W and D denotes molten pool width and molten pool depth, respectively.
As shown in Figure 8, the closest agreement was observed with 6 × 6 heat sinks under test 1 conditions. The horizontal axis 0 × 0 denotes the calculated molten pool dimensions without a heat sink. In other words, the heat loss from the part boundary was not considered. Expt. denotes the experimental molten pool dimensions. The more heat sinks, the smaller the molten pool dimensions, which is consistent with the heat-affected zone dimensions. This trend confirms the instinctive trend that the more heat loss, the smaller the heat-affected zone and molten pool, and vice versa. which is consistent with the heat-affected zone dimensions. This trend confirms the instinctive trend that the more heat loss, the smaller the heat-affected zone and molten pool, and vice versa. The calibration data from post-process temperature measurement is more reliable and less sensitive to operation than the in-situ measurements from thermocouples and infrared cameras [10,11]. The experimental molten pool measurements were conducted at least in triplicate with negligible deviation observed under each process condition. Therefore, the use of post-process measurement of molten pool dimensions for calibration purposes is acceptable. The experimental calibration process is necessary for the presented model to avoid the compensation of computational The calibration data from post-process temperature measurement is more reliable and less sensitive to operation than the in-situ measurements from thermocouples and infrared cameras [10,11]. The experimental molten pool measurements were conducted at least in triplicate with negligible deviation observed under each process condition. Therefore, the use of post-process measurement of molten pool dimensions for calibration purposes is acceptable. The experimental calibration process is necessary for the presented model to avoid the compensation of computational efficiency using iterative calculation-based calibration on heat sink temperatures. The optimal heat sink settings are material-dependent because the material-dependent heat transfer coefficients of convection and radiation are used in the presented model.
Moreover, the temperature calculations were carried out with a MATLAB (MathWorks, USA) program on a personal computer running at 2.8 GHz. The part size was 1 mm × 1 mm × 0.2 mm, with increments of 2.5 µm in all directions. The computational times with 4 × 4, 6 × 6, 8 × 8, and 10 × 10 heat sinks were 49.42 s. 94.44 s, 160.28 s, and 239.62 s, respectively. For comparison, FEM models required at over an hour for similar calculation increments and prediction accuracies [12,13,16,17]. The semi-analytical models [38,39] using coarse mesh resolution for heat transfer boundary conditions will significantly reduce the prediction accuracy for large-scale parts. The presented model does not rely on FEM or any iteration-based calculation, and thus remains unaffected.
With 6 × 6 heat sink settings, the molten pool dimensions were predicted under 8 different process conditions, as shown in Table A2 in the Appendix A. The predicted molten pool dimensions were validated by the experimental measurements. The continuity of the scan tracks was confirmed from the observation on the top view. The experimental measurement was made at least in triplicate with negligible variation observed under each process condition. Close agreement was observed, as shown in Figure 9. The deviation might be caused by the simplified point heat source solution without considering the heat source profile, the adopted absorptivity, and material properties. The associated data are given in Table A2 in the Appendix A. With 6 × 6 heat sink settings, the molten pool dimensions were predicted under 8 different process conditions, as shown in Table A2 in the Appendix A. The predicted molten pool dimensions were validated by the experimental measurements. The continuity of the scan tracks was confirmed from the observation on the top view. The experimental measurement was made at least in triplicate with negligible variation observed under each process condition. Close agreement was observed, as shown in Figure 9. The deviation might be caused by the simplified point heat source solution without considering the heat source profile, the adopted absorptivity, and material properties. The associated data are given in Table A2 in the Appendix A. The presented model has demonstrated high prediction accuracy and high computational efficiency using the closed-form solutions without resorting to FEM or any iteration-based calculations. Those advantages allow the presented model to be used for temperature prediction for large-scale parts and process-parameter planning through inverse analysis [45,46]. In the future, the presented model should be further developed for temperature prediction of multi-track and multilayers scans. The applicability of the developed model on other widely used powder materials should The presented model has demonstrated high prediction accuracy and high computational efficiency using the closed-form solutions without resorting to FEM or any iteration-based calculations. Those advantages allow the presented model to be used for temperature prediction for large-scale parts and process-parameter planning through inverse analysis [45,46]. In the future, the presented model should be further developed for temperature prediction of multi-track and multi-layers scans. The applicability of the developed model on other widely used powder materials should be investigated. The intensity profile of the laser heat source should also be considered as a further improvement. The commonly observed defects of residual stress, porosity, and part deviation can be further investigated based on the calculated temperatures and temperature gradient because those defects were caused by the elevated temperature levels and large temperature gradients.

Conclusions
This work presented an analytical model for temperature prediction in powder bed metal additive manufacturing (PBMAM), also known as powder bed fusion (PBF) or selective laser melting (SLM). This model was developed based on physics with consideration of heat conduction, convection, and radiation, heat absorption, and latent heat affect. It was constructed from two closed-form solutions, namely the moving point heat source solution and the heat sink solution. The original heat sink solution was developed based on the heat transfer equations for convection and radiation. The influence of the chosen number of heat sinks on the predicted temperature profiles, temperature spatial gradient, and temperature-affected material property variation were investigated. An optimal number of heat sinks was determined by calibration with documented experimental values. The presented model was validated against documented experimental values under different process conditions. Close agreements were observed upon validation.
The presented model has promising short computational time without resorting to FEM or any iteration-based simulations, which was confirmed from the recorded computational time. With the heat sink solution, the complete understanding of the heat transfer mechanism in PBMAM can be implemented effectively and efficiently. The employed heat sink solution improves the prediction accuracy, and thus the usefulness of the analytical modeling in real applications, specifically the temperature investigation in PBMAM. It can be used for temperature prediction of large-scale parts and process-parameter planning through inverse analysis because of the high computational efficiency. The calculated temperature and material property variation allows further investigations of residual stress, porosity, and part deviations, which are caused by repeated rapid heating and solidification.
Author Contributions: J.N. performed the formal analysis and investigation, extracted, analyzed, and validated the data, and wrote the manuscript. S.Y.L. supervised the research project and proofed the manuscript. D.E.S. and H.G. provided resources and critical feedback.
Funding: This research project was funded by the Boeing company. The grant number is confidential.