## Abstract

Temperature history prediction is essential for a better understanding of the relationship between microstructural change and processing conditions for energy beam additive manufacturing fabricated components. Here, a new efficient approach combining a moving heat source analytical model with a melting and solidification model is presented. An innovative method is proposed to compute the “effective computation zone” as a boundary condition, which can save computation time significantly. Notably, the computational efficiency can improve by 10^{4}–10^{5} compared with finite element models. With this range of improvement efficiency, the temperature predicted based on this method is consistent (around 9% of average deviation) with experimental measurements by the thermocouple. This model can be used as a reference to define the boundary condition for further complex numerical analysis with improved accuracy at a reduction of efficiency as desired. In addition, it can be used as a reference to determine processing conditions that would allow the efficient and effective control of the temperature history within a range for a certain microstructure design.

## 1 Introduction

Additive manufacturing has its inherent cyclic heating and cooling process that would affect the microstructure and properties of components fabricated. As a result, it is crucial to use heat transfer models to compute the temperature history for processing control. Based on a thermal model, many researchers have investigated phase transformation, microstructure variation, defects and porosity, and changes in mechanical properties due to the different processing conditions of energy beam additive manufacturing [1–8].

Modeling the temperature history of the additive manufacturing process is complex and time-consuming. Lee and Zhang presented a discrete element model (DEM) to compute the temperature for a single layer over a substrate of 1 mm × 0.4 mm; it took less than 40 h of CPU clock time to finish the simulation [9]. Despite the small domain, the time cost of the finite element model (FEM) is high. Fachinotti and Cardona proposed a three-dimensional (3D) FEM for a thin wall (3 mm × 38 mm × 13 mm) with a CPU computation time of around 5288 s [10]. Panagiotis presented a transient 3D FEM for a domain of 20 cm × 10 cm × 10 cm, which required 6296 s for computation [11]. Other models such as the level set, the volume of fluid with a finite difference, lattice Boltzmann, and arbitrary Lagrangian–Eulerian methods are even more computationally intensive. Morville et al. reported that it took around 4–5 days to simulate the temperature for a single layer (35 mm long) and around 4 weeks for five layers [12]. Khairallah and Anderson presented a 3D mesoscopic model for a solid substrate (1000 *µ*m × 300 *µ*m × 50 *µ*m), which required 100,000 h of CPU computation time [13]. Therefore, it is computationally expensive to obtain an accurate prediction of temperature history. Yang et al. developed a semi-analytical model that can combine the analytical solution with fine-meshed FEM as a boundary condition. However, it is still computationally expensive since the model took 7–15 min for the case of only a single layer. Ning et al. developed an analytical model for powder bed additive manufacturing. However, they did not propose any potential methods to save computational time, which is one of the key aspects of using an analytical model [14]. The purpose of this paper is to fill the gap of computational cost by using the concept of “effective computational zone.” An efficient heat transfer model would be an ideal tool for the rough estimation of the temperature with an acceptable compromise in prediction accuracy, which is useful especially in the early stage of designing manufacturing processing conditions.

The analytical model has advantages including low-cost computational, simplicity, and ease of use. Various studies have proposed a moving heat source model that can be used for energy beam additive manufacturing. Carslaw and Jaeger presented a method to solve a stationary heat source problem mathematically [15]. Rosenthal proposed the solution by integrating the temperature response of an infinite number of subsequent instantaneous points [16]. Hou and Komanduri developed a model that can be used for both transient and steady-state conditions. In addition, a general solution was provided for different shapes and different intensity distributions of the heat source [17]. Eagar and Tsai introduced a parameter as a standard deviation for a Gaussian distribution two-dimensional (2D) heat source [18]. Goldak et al. proposed a double-ellipsoidal method to separate the heat source into the rear and front end [19]. Furthermore, Nguyan et al. developed a closed-form analytical model to simulate a 3D moving heat source [20]. Yang et al. developed a semi-analytical model that computes the heating area by an analytical model and corrects the boundary field by finite difference [21]. Li et al. considered the case of multi-layers directed energy deposition by virtual heat source method [22]. However, a central issue for these models is they did not consider the melting and solidification process, which is critical to the additive manufacturing processes that involve melting of powders or raw materials such as the selective laser melting process. The latent heat required for the phase change between solid and liquid phase can affect the temperature computed, especially for additive manufacturing that involves multi-heating and multi-cooling processes.

The objective of this study is to develop an approach that can combine an analytical model with melting/solidification models to compute the temperature history efficiently with acceptable error for energy beam additive manufacturing.

## 2 Modeling Approaches

A problem of energy beam additive manufacturing the interactions between multiple layers with the melting and solidification processes and a moving heat source. An analytical model is presented to compute the temperature rise and fall due to a relative material position to the heat source. Figure 1 shows the flowchart of the analytical model. First, the effective computation zone (ECZ) was defined based on the accuracy needed. Second, the heat sources within the effective computation zone both were recorded currently and previously. Third, the surface temperature for the point of interest was computed based on the moving heat source model. Fourth, if the surface temperature was equal to or higher than the melting point, the melting and solidification process was considered as a moving boundary problem with a modified model that includes latent heat and phase transformation. Fifth, the temperature changes due to the previous heat sources were computed along with the cooling process. Sixth, the computed temperature is the combination of temperature change due to both the current heat source and the previous heat source based on superposition. Seventh, the physical properties and solidified layers for each time-step were updated for the next computation loop.

### 2.1 Modeling of the Temperature Distribution Due to a Moving Heat Source.

*T*is the temperature;

*x*,

*y*, and

*z*are coordinates along with the three mutually orthogonal directions;

*α*is the thermal diffusivity,

*t*is the time,

*ρ*is the density,

*C*is the heat capacity coefficient, and

_{p}*q*is the heat source.

*P*is the total energy released from the heat source,

*C*

_{p}is the heat capacity,

*ρ*is the density of the medium, and

*r*is the distance between the relevant position and the heat source, which is

*r*= [(

*x*−

*x*

_{i})

^{2}+ (

*y*−

*y*

_{i})

^{2}+ (

*z*−

*z*

_{i})

^{2}]

^{1/2}.

*E*,

*F*,

*G*, and area

*A*

_{pl}are listed in Table 1.

Shape of heat source | Distribution of intensity | E | $F$ | $G$ | $Apl$ |
---|---|---|---|---|---|

Elliptical (or circular disc) | Uniform | 1 | 1 | 1 | πa_{0}b_{0} (or $\pi r02$) |

Parabolic | 0.5 | $(1\u2212n2)$ | $(1\u2212m21\u2212n2)$ | ||

Normal | 0.1079 | exp[−(3n)^{2}] | $exp[\u2212(3m1\u2212n2)2]$ | ||

Rectangular (or square) | Uniform | 1 | 1 | 1 | 4a_{0}b_{0} (or 4$a02$) |

Parabolic | 4/9 | 1 − n^{2} | $1\u2212m2$ | ||

Normal | π/36 | $exp[\u2212(3n)2]$ | $exp[\u2212(3m)2]$ |

Shape of heat source | Distribution of intensity | E | $F$ | $G$ | $Apl$ |
---|---|---|---|---|---|

Elliptical (or circular disc) | Uniform | 1 | 1 | 1 | πa_{0}b_{0} (or $\pi r02$) |

Parabolic | 0.5 | $(1\u2212n2)$ | $(1\u2212m21\u2212n2)$ | ||

Normal | 0.1079 | exp[−(3n)^{2}] | $exp[\u2212(3m1\u2212n2)2]$ | ||

Rectangular (or square) | Uniform | 1 | 1 | 1 | 4a_{0}b_{0} (or 4$a02$) |

Parabolic | 4/9 | 1 − n^{2} | $1\u2212m2$ | ||

Normal | π/36 | $exp[\u2212(3n)2]$ | $exp[\u2212(3m)2]$ |

*x*-direction as an example, the temperature can be modeled in Eqs. (4)–(6). The heat source has a length and width of 2

*a*

_{0}and is moving constantly along the

*x*-direction as shown in Fig. 2.

*T*

_{hs}is the surface temperature of point (

*x*,

*y*,

*z*) interest, (

*x*

_{i},

*y*

_{i},

*z*

_{i}) is the moving heat source position as a function of time,

*v*is the scanning speed,

*α*is the thermal diffusivity, and

*P*is the power of heat source.

### 2.2 Modeling of the Melting and Solidification Processes.

Given the surface temperature computed with Eq. (3), as the temperature of the surface reaches the melting point, the melting and solidification process can be simulated as a moving boundary problem. The boundary in this instance is the interface of the solid phase and liquid phase. The moving boundary problem has been investigated in previous studies, and a solution was proposed by Stefan [23]. In this study, the model assumed that the thermal properties were constant for a semi-infinite slab and that it was heated at a constant temperature that is higher than the melting point of one end. The interface between the solid and liquid was a moving boundary designated as *S*(*t*) as shown in Fig. 3. In our study, we revise Stefan's moving boundary problem with the heat source on the top surface as shown in Fig. 4. The phase change is assumed to occur at a single temperature without the mixture of a two-phase zone.

The transient moving boundary problem can be solved mathematically by a heat diffusion equation with initial and boundary conditions in Eqs. (7)–(9):

*H*is the latent heat,

*T*

_{m}is the melting temperature,

*T*

_{hs}is the surface temperature, and

*k*is the thermal conductivity. Neumann [24] proposed an analytical solution for Eqs. (7)–(9). The solution was based on the error function as shown in Eq. (10):

*α*

_{l}is the thermal diffusivity for the liquid phase,

*α*

_{s}is the thermal diffusivity for the solid phase,

*erf*is the error function, and

*φ*is the solution for Eqs. (11) and (12) [23,25,26]:

*Ste*

_{l}is the Stefan number for the liquid phase and

*Ste*

_{s}is the Stefan number for the solid phase.

### 2.3 Defining Effective Computation Zone.

The ECZ is the defined area in which the material points and heat source within the ECZ would affect the temperature of the point concerned. Each material point is the point in the physical modeling/simulation domain. The ECZ cutoff is based on the tolerance temperature defined. The material points outside of the effective computation zone do not need to be considered when computing the point concerned. As the acceptable error decreases, the effective computation zone and time required for computation increase, but the accuracy improve. When applying the concept of an effective computation zone, the boundary of the effective computation zone is set to be room temperature. In addition, when computing temperature cooling from the previous time, only the heat source within the effective computation zone needs to be considered. This will be further discussed in Sec. 2.4.

### 2.4 Modeling Temperature Change Due to Previous Heat Source.

*z*= 0) can transfer heat efficiently so that the temperature maintains as room temperature constantly as shown in Fig. 5.

*x*, Δ

*y*, and Δ

*z*are the distance between the point concerned and boundary of effective computation zone shown in Fig. 5,

*A*is the contact area of each element,

*V*is the volume of the element, Δ

*t*is the cooling time,

*T*

_{c}is the temperature before the cooling process, and

*T*

_{c, new}is the temperature after the cooling process.

*T*

_{h}) and previous moving heat source (

*T*

_{c}) written in Eq. (17) based on superposition. This algorithm is shown in Fig. 1.

## 3 Results and Discussion

The presented model in Sec. 2 was coded in matlab software to compute the temperature due to the moving energy beam. In order to validate the model, the computational result was compared with the experimental measurement and FEM results (Sec. 3.1). The computational efficiency was compared with FEM and DEM (Sec. 3.2), which demonstrated the power of the presented analytical model.

### 3.1 Model Validation.

An experiment was conducted by Amine et al. [27] focusing on a thin, multilayer, 316L stainless steel wall fabricated by a 600 W laser. The thin wall was composed of 87 layers (length of 30 mm, total height of 14 mm, and width of 2.5 mm), as shown in Fig. 6. The laser diameter was 5 mm, and the scanning strategy was zigzag at a constant travel velocity of 450 mm/min. The cooling substrate was 50 mm wide, 50 mm long, and 12.7 mm thick. In order to measure the temperature, a K-type thermocouple was setup at the end of the laser track beneath the substrate at a 3 mm depth. The chemical composition of 316L stainless steel is as shown in Table 2. The thermal conductivity used in the simulation was temperature-dependent and was updated at each time-step as shown in Table 3. The constant material properties are listed in Table 4. The model assumed that the raw materials were fully melted and had no porosity.

Element | C | Mn | P | Si | Ni | Cr | Mo | Fe |
---|---|---|---|---|---|---|---|---|

(W%) | <0.03 | <2 | <0.045 | <1 | 10–14 | 16–18.5 | 2–3 | Bal |

Element | C | Mn | P | Si | Ni | Cr | Mo | Fe |
---|---|---|---|---|---|---|---|---|

(W%) | <0.03 | <2 | <0.045 | <1 | 10–14 | 16–18.5 | 2–3 | Bal |

Temperature, K | 300 | 400 | 500 | 600 | 700 | 800 | 1000 |
---|---|---|---|---|---|---|---|

Thermal conductivity, W/m K | 13.4 | 15.2 | 16.75 | 18.3 | 19.8 | 21.3 | 24.2 |

Temperature, K | 300 | 400 | 500 | 600 | 700 | 800 | 1000 |
---|---|---|---|---|---|---|---|

Thermal conductivity, W/m K | 13.4 | 15.2 | 16.75 | 18.3 | 19.8 | 21.3 | 24.2 |

Melting temperature, K | 1385 |

Latent heat, kJ/kg | 272.5 |

Density, kg/m^{3} | 8000 |

Specific heat for solid, $J/kgK$ | 500 |

Specific heat for liquid, $J/kgK$ | 820 |

Boiling temperature, K | 3900 |

Laser absorptivity | 0.6 |

Melting temperature, K | 1385 |

Latent heat, kJ/kg | 272.5 |

Density, kg/m^{3} | 8000 |

Specific heat for solid, $J/kgK$ | 500 |

Specific heat for liquid, $J/kgK$ | 820 |

Boiling temperature, K | 3900 |

Laser absorptivity | 0.6 |

To further validate the accuracy of the analytical model, another 3D transient FEM is presented using abaqus/cae software [27]. The meshes were built using hypermesh software with eight-node linear brick-type elements. Each element was around 0.5 mm × 0.5 mm × 0.5 mm.

Figure 7 shows temperature measurements by a thermocouple compared with simulation results (including the analytical model and FEM results). Since the thermocouple was 3 mm below the substrate, the temperature was lower than the melting point to ensure the accuracy of the thermocouple. The zigzag pattern of the increasing and decreasing temperature computed may be attributed to the relative position between the moving heat source and the position during the scanning process. The variation in heating and cooling was larger at the beginning of the fabrication process, and the temperature reached a steadily increasing state as the layers were built. This change may be attributed to the cooling process by a substrate. As the height of building layers increased, there would be increased thermal resistance to the transfer of heat to the bottom of the cooling substrate, leading to the storage of heat within the thin wall.

*ɛ*

_{1}). On the other hand, if we increase the ECZ like model 3, the error term (

*ɛ*

_{3}) can be decreased. Figure 9 illustrates that the accuracy can be improved by increasing the effective computation zone and is proportional to the time required for computation as shown in Table 6. The accuracy can be further improved by changing the effective computation zone at a reduction of efficiency as desired. The temperature computed by the presented analytical model was higher than that obtained by experimental measurement, which may be associated with the assumption of a semi-infinite cooling substrate and the cutoff of an effective computation zone. The accuracy can be further improved by changing the effective computation zone at a reduction of efficiency as desired.

Model | Average deviation | Effective computation zone radius |
---|---|---|

FEM | 23% (underestimated) | |

Proposed model 1 | 24% (overestimated) | 0.5 mm |

Proposed model 2 | 13% (overestimated) | 1 mm |

Proposed model 3 | 8% (overestimated) | 2 mm |

Model | Average deviation | Effective computation zone radius |
---|---|---|

FEM | 23% (underestimated) | |

Proposed model 1 | 24% (overestimated) | 0.5 mm |

Proposed model 2 | 13% (overestimated) | 1 mm |

Proposed model 3 | 8% (overestimated) | 2 mm |

Model | Effective computation zone radius | Required CPU computation time |
---|---|---|

Proposed model 1 | 0.5 mm | 7.582 s |

Proposed model 2 | 1 mm | 11.321 s |

Proposed model 3 | 2 mm | 15.527 s |

Model | Effective computation zone radius | Required CPU computation time |
---|---|---|

Proposed model 1 | 0.5 mm | 7.582 s |

Proposed model 2 | 1 mm | 11.321 s |

Proposed model 3 | 2 mm | 15.527 s |

### 3.2 Melt Pool Distribution.

The melting and solidification process was computed based on the moving boundary of the solid–liquid interface of the Neumann solution [24]. Figure 10 shows the melt pool distribution, in which the melting boundary is penetrated from the surface to about 21.6 *µ*m in depth. The simulation result shown in Table 7 is compared with the finite element model [8], which was already validated by the experiment. The predicted melt pool dimension based on the presented model was 68.0 *µ*m in width (about 8% difference) and 21.6 *µ*m in width (about 14% difference) compared with the finite element model. The deviation is defined in Eq. (18), and the reference is the melt pool predicted by FEA.

### 3.2 Computational Efficiency and Model Application.

The presented analytical model was computationally less intensive, simplified, and easy to use. The computational efficiency can be significantly improved compared with that of the FEM and DEM. Figure 11 shows a comparison of the CPU time required for computation. As the computed elements were increased for the FEM and DEM, the computational time was increased to recombine all sets of element equations into the entire boundary condition. On the other hand, the analytical solution computed only the element within the effective computation zone based on the relative position of the relevant point to the heat source.

Table 8 and Fig. 12 show a comparison of the computational time required for a thin plate of 1 m × 1 m × 0.25 mm (10,000,000 elements). Although there were different processing conditions, CPUs, and numbers of cores used in each case, the comparison gives a rough estimation of how efficiency improves with our analytical model. In Table 8, the efficiency can be compared because the computational efficiency was based on the average CPU time per element and per timestep computed. It is indicated that the proposed model could improve efficiency by 10^{5} to 10^{6} for 10,000,000 elements.

CPU time, s | Proposed model efficiency improvement | CPU | # of cores | Reference | |
---|---|---|---|---|---|

Proposed model | 2.26 × 10^{1} | Intel^{®} Core^{™} i7-8550U CPU, 1.8 GHz 2.00 GHz | 4 | ||

2D FEM | 4.89 × 10^{5} | 2.16 × 10^{4} | Intel Core 2 Duo 6300 | 2 | [10] |

DEM | 6.91 × 10^{5} | 3.06 × 10^{4} | Intel^{®} Xeon^{®} Processor E5335 4 GB RAM | 4 | [9] |

3D FEM | 7.02 × 10^{5} | 3.10 × 10^{4} | 3.1GHz | 16 | [11] |

3D FEM | 3.90 × 10^{6} | 1.72 × 10^{5} | Intel Core 2 Duo 6300 | 2 | [10] |

CPU time, s | Proposed model efficiency improvement | CPU | # of cores | Reference | |
---|---|---|---|---|---|

Proposed model | 2.26 × 10^{1} | Intel^{®} Core^{™} i7-8550U CPU, 1.8 GHz 2.00 GHz | 4 | ||

2D FEM | 4.89 × 10^{5} | 2.16 × 10^{4} | Intel Core 2 Duo 6300 | 2 | [10] |

DEM | 6.91 × 10^{5} | 3.06 × 10^{4} | Intel^{®} Xeon^{®} Processor E5335 4 GB RAM | 4 | [9] |

3D FEM | 7.02 × 10^{5} | 3.10 × 10^{4} | 3.1GHz | 16 | [11] |

3D FEM | 3.90 × 10^{6} | 1.72 × 10^{5} | Intel Core 2 Duo 6300 | 2 | [10] |

The results revealed that the computational efficiency was significantly improved. The large computation time difference may be attributed to the cutoff temperature and time needed for computation. In the FEM, the whole processing region needs to be computed based on the boundary condition to obtain the temperature for a specific region. On the other hand, the analytical model focused on the relevant region and predicted the temperature based on the effective computation zone and time.

The proposed closed-form analytical model provides a tool that can efficiently compute the temperature history because it deploys a moving heat source during the energy beam additive manufacturing process. Figure 12 shows the model application: the processing conditions of manufacturing such as scanning velocity can be related to material structural change. Based on the temperature history, the critical temperature that leads to phase transformation can be determined. Furthermore, based on the cooling rates from the liquidus to solidus temperature, the microstructure feature and material properties after solidification can be further investigated [1–8]. The strategy of using an effective computation zone can reduce the computation cost significantly and further save the design iteration cycle time.

For future studies, the proposed model can be generalized by substitution with different moving heat sources such as square, uniform, and point heat sources. The boundary condition, dimensions, temperature, and heat flux levels that can be treated are based on the analytical model applied. For example, if we used Cheng and Lin [28] model, the heat flux is Gaussian distribution but is a circle 2D heat flux. If we applied the Nguyan analytical model [20], we can simulate a 3D moving heat source condition. This paper illustrated a method to combine the melting/solidification model with the analytical heat source model to get an accurate computational result. We can also combine more sophisticated melting/solidification models with the analytical heat source model selected by the same approach stated in this paper.

The model assumed that the material properties are linear systems so that the superposition method can be applied. However, to handle nonlinearity in most practical cases, we can assume that the property at each discrete time-step is linear and integrate the performance of each time-step to simulate the overall nonlinear results. The different substrate geometries, layers geometries, and material properties can be considered by integrating and updating the parameters at each time step to improve the accuracy. For example, we can include the heat transfer of solidified printed layers by updating the distance between the surface heat source and point concerned as shown in Fig. 13. Another example is the nonlinear thermal conductivity of 316L stainless steel at a different temperature as demonstrated in Sec. 3.1. The thermal conductivity used in the simulation was temperature-dependent and was updated at each time-step to include the nonlinearity. In addition, the phase change temperature can be modeled as a range between the liquidus temperature [29] and solidus temperature to improve the accuracy by using the same method.

## 4 Summary and Conclusions

In summary, an efficient approach is presented to compute the temperature history of components fabricated by energy beam additive manufacturing. The method combined an analytical model with a melting/solidification model. The predicted temperature distribution is consistent (around 9% average deviation) with thermocouple measurements, and the accuracy was better than the FEM simulation [27]. The computational efficiency is improved 10^{4}–10^{5} compared with that of FEM and DEM for a thin plate case study of a 1 m × 1 m × 0.25 mm (10,000,000 elements) thin plate. The considerably shorter computation time has attributed to the cutoff of the area required for defining the boundary condition needed for computation based on the concept of an effective computation zone. The presented model provides an efficient tool for temperature computation that is less computationally intensive for energy beam additive manufacturing.

## Acknowledgment

This work was supported by the National Science Foundation (CMMI), through Grant No. 1562960. Some ideas and materials are patent pending owned by Purdue University.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper.