Fast Detection of Heat Accumulation in Powder Bed Fusion Using Computationally Efficient Thermal Models

The powder bed fusion (PBF) process is a type of Additive Manufacturing (AM) technique which enables fabrication of highly complex geometries with unprecedented design freedom. However, PBF still suffers from manufacturing constraints which, if overlooked, can cause various types of defects in the final part. One such constraint is the local accumulation of heat which leads to surface defects such as melt ball and dross formation. Moreover, slow cooling rates due to local heat accumulation can adversely affect resulting microstructures. In this paper, first a layer-by-layer PBF thermal process model, well established in the literature, is used to predict zones of local heat accumulation in a given part geometry. However, due to the transient nature of the analysis and the continuously growing domain size, the associated computational cost is high which prohibits part-scale applications. Therefore, to reduce the overall computational burden, various simplifications and their associated effects on the accuracy of detecting overheating are analyzed. In this context, three novel physics-based simplifications are introduced motivated by the analytical solution of the one-dimensional heat equation. It is shown that these novel simplifications provide unprecedented computational benefits while still allowing correct prediction of the zones of heat accumulation. The most far-reaching simplification uses the steady-state thermal response of the part for predicting its heat accumulation behavior with a speedup of 600 times as compared to a conventional analysis. The proposed simplified thermal models are capable of fast detection of problematic part features. This allows for quick design evaluations and opens up the possibility of integrating simplified models with design optimization algorithms.


Introduction
Additive manufacturing (AM) offers unprecedented design freedom as compared to conventional manufacturing techniques. The layer-by-layer material deposition allows for manufacturing functional parts with high geometric complexity. Due to this advantage, AM has already gained significant popularity among manufacturing industries such as automotive, aerospace and medical [1,2]. However, the AM process is not free from manufacturing constraints and, if overlooked, these constraints can cause a wide range of defects in the final part which lead to an increased overall cost. Therefore, to fully capitalize on the benefits offered by AM, manufacturing limitations should be considered at the design stage.
Laser powder bed fusion (LPBF) process is one of the most common techniques for printing metal parts. It involves sequential deposition of metal powder layers which are selectively molten and fused together in predefined areas using a moving laser beam. This implies that heat primarily flows from each newly deposited topmost layer towards the thick baseplate at the bottom, which acts as a heat sink. It is observed that whenever the deposited heat is not properly evacuated towards the baseplate, it leads to local overheating or heat accumulation [3]. It typically refers to a situation where material locally experiences thermal process conditions that result in temperatures outside the desired temperature range needed to obtain the desired final product quality. Given the wide range of potential temperature-induced production failure mechanisms, overheating can manifest itself in many forms. In the in situ monitoring studies conducted by Hooper [4] and Kruth [5], overheating is characterized by the enlarged melt pool observed in the vicinity of lower conductivity regions which obstruct heat flow. In the numerical study conducted by Hodge et al. [6], the overheating phenomenon is characterized by the overshoot of simulation temperatures above the melting point.
It has been widely recognized in the literature that local overheating or heat accumulation is detrimental for final part quality and can cause various defects. For example, Mertens et al. [3] investigated the defect of dross formation where an enlarged melt pool caused by local overheating leads to undesired sintering of loose powder in its vicinity. Charles et al. [7] also studied dross formation by investigating the correlation between process parameters, e.g., laser power, scan speed and the resulting surface roughness. High process temperatures along with the effect of surface tension lead to the defect of melt ball formation [8]. These cases are examples where overheating depends on the local distribution of conductivity in the vicinity of melt pool. Another more recently studied phenomenon is that of gradual accumulation of heat with increasing part height. As more layers are deposited, they act as a thermal barrier causing prolonged heated zones with slower cooling. This becomes increasingly significant in tall parts with high heat capacity. Therefore, the phenomenon was initially investigated for processes such as laser metal deposition (LMD) [9] and wire-arc-based AM process (WAAM) [10], where typical part sizes are larger than LPBF. More recently, it has been studied for the LPBF process as well. For example Hilgenberg [11] used in situ thermography for observing cooling conditions during the process. It is reported that heat accumulation increases with increasing build height. Also, inter-layer time (ILT), i.e., the time elapsed between processing successive layers, is found to be an important parameter which inversely affects heat accumulation.
It is important to understand that both these phenomena, i.e., local overheating and gradual heat accumulation during the build, are intimately coupled. Elevated temperatures due to local overheating contribute to the gradual accumulation in the successive layers and accumulated heat will conversely intensify local overheating. Therefore, both these phenomena jointly determine the overall thermal history of the part, which subsequently has a direct effect on the final part quality. For example, Hilgenberg [11] performed metallographic analysis to study the influence of heat accumulation on observed microstructures. It is suggested in this paper that specific temperature ranges which are important from the context of phase transformations should be considered while analyzing the influence of heat accumulation. Another aspect investigated by Wildman [12] shows that regions susceptible to overheating generate significant residual stresses, which are likely to result in undesired deformations and potentially build failures. Kastsian and Reznik [13] analyzed a part geometry where local overheating led to excessive deformations causing recoater jamming. Heigel and Whitenton [14] presented layer-wise cooling rate maps which show that overheating zones tend to cool slower. Lastly, as suggested by Charles et al. [7], another aspect of overheating control is associated with selecting optimal process parameters. For example, Mertens et al. [3] varied laser power in accordance with the level of overhang for controlling dross formation. However, such approaches also rely on prior knowledge of overheating zones so that parameters can be tuned accordingly. Hence, early identification of geometric features which are prone to local overheating can assist designers and machine operators in judging manufacturability of the final design.
Typically, downfacing or overhanging features are prone to this phenomenon as loose powder, with significantly lower thermal conductivity than that of the bulk material, does not allow for rapid heat evacuation. The degree of overheating typically increases with decreasing the overhang angle defined as the angle between part surface and the base plate. Therefore, heuristic AM design guidelines are used, where features with overhang angles smaller than a critical value are avoided [15,16]. Although overhangs are a salient example of features which cause overheating, the phenomenon is not uniquely linked to them. For example, LPBF specimens manufactured by Zimm. [17] and Patel et al. [18] demonstrate overheating induced discoloration and high surface roughness, respectively, even after strictly following the overhang design guideline. This indicates that not only the local overhang angle, but the thermal response should be considered for detecting and preventing overheating.
Numerical modeling of the transient heat transfer phenomena during the LPBF process can provide thermal history information which can be used for identification of features causing local overheating. However, excessive computational cost associated with detailed modeling of the LPBF process prohibits part-scale implementation, particularly in iterative design settings. Please note that the size of the laser spot is in the order of a few microns whereas the part size is in the order of centimeters. Moreover, in time domain, the phenomenon of melting and re-solidification is extremely fast, whereas cooling between the layers is in order of tens of seconds. Hence, to capture the steep thermal gradients caused by a fast-moving laser, very fine spatial and temporal resolution is required which makes the modeling process computationally intractable. To address this issue, researchers reduce the computational complexity of models by introducing simplifications in various ways. For example, several studies motivate consideration of conduction as the only mode of heat transfer while neglecting convection and radiation [19][20][21]. Another commonly used simplification is to simulate an entire layer addition as a single simulation step, i.e., the entire layer area is simultaneously exposed to a heat flux. It has been demonstrated by Rubenchik [22] that during a new layer deposition, 4-8 previously deposited layers are also re-molten and fused. This re-melting is intentional as it allows for a seamless connection between the layers. Therefore, the concept of "superlayers" or layer lumping has been proposed where deposition of one superlayer refers to the simultaneous deposition of multiple real layers while modeling the process [23][24][25][26].
These simplifications can provide significant computational gains as compared to high-fidelity models, with a reasonable compromise in accuracy. For example, Davies [24] studied the impact of excluding powder conduction, simplifying laser scan strategies and lumping multiple layers during thermal simulation. The simulation data is analyzed for accuracy by comparing it with the experimental temperatures recorded using thermocouples embedded inside the printed parts. It is reported that powder exclusion results in 4 times faster simulation. Furthermore, the layer-by-layer approach is reported to be 6 times faster than the model which simulates scan strategies. Lastly, lumping of 4 layers provided another 3 fold decrease in the computational time. It is reported that with increased layer lumping, the model captures an average evolution of the temperature field during the manufacturing process. Another study by Harrison [26] investigated the effect of using lumped layers as thick as 4, 12 and 24 times the real layer thickness. The analysis was carried out for a 2D real-size part (height = 54.7 mm, width = 20.4 mm) and it was reported that lumping of 4 layers already results in a 20 times faster simulation compared to the one with no layer lumping. In both of these cases, although approximations provide significant computational gain, still the simulation times remain in the order of hours. For example, out of the two cases discussed above, the former reports that simulation with 4 layer lumping requires 1.33 h for parts with volume of 107 cm 3 printed with 992 layers. For the latter 2D example, the model with lumping of 4 layers needed 4.6 h to complete.
In this paper, we aim to develop a model which could identify zones of local overheating and in order to assess manufacturability of designs and make quick design decisions, even faster computation times will be preferred. This is also motivated by the desire of integrating AM models within an optimization framework. Therefore, in this research, first a lumped layer-by-layer model is considered for simulating a real-size LPBF part and the temperature data is used to judge the heat accumulation tendencies of different design features. Details about this model are presented in Section 2. As highlighted in the previous paragraph, even with the simplifications of powder exclusion and lumped layer-by-layer deposition technique, the simulation time remains prohibitive for large scale implementations. Therefore in Section 3, six additional simplifications are presented. The first three have been partially discussed in the literature but an analysis of their impact on specifically overheating detection was never done before, which is discussed in Section 3.1. Next, three further simplifications are proposed in Section 3.2 which are the primary novel contributions of this paper. An analysis based on the analytical solution of the one-dimensional heat equation provides the motivation for these novel simplifications. Subsequently, Section 4 presents results obtained using the different simplifications and it is shown that the proposed simplifications provide adequate prediction of overheating regions while providing significant computation gains. Lastly, conclusions and future directions are discussed in Section 5.

Reference Thermal LPBF Process Model
In this section, a layer-by-layer part-scale thermal model is presented, intended for judging manufacturability of designs from the context of local overheating. It is based on the thermal model previously presented by Chiumenti et al. [27] for simulating material deposition in the laser material deposition (LMD) process. The same concept was later used for LPBF modeling in several publications [24,28,29]. The experimental validation of the model has been presented in Davies [24] where simulation results were compared with temperatures empirically recorded using in situ thermocouples placed inside the part during the LPBF process. It was shown that this modeling method can correctly capture the thermal history as obtained using the thermocouples. More importantly, it was observed by Neiva et al. [29] that the model correctly predicts the high temperatures measured near an overhanging feature. Therefore, it is deemed suitable for the intended purpose of identifying geometric features which cause local heat accumulation during the LPBF process. The model presented in this section acts as a reference for the simplified models presented in the later sections.

Model Description
A typical LPBF process includes layer-by-layer deposition of material on a thick baseplate. Figure 1 shows a schematic where an arbitrary 3D object is considered with volume Ω already deposited. The surfaces ∂Ω of the partly manufactured object are classified as top, lateral and bottom represented as ∂Ω top , ∂Ω lat and ∂Ω bot , respectively. The part remains submerged inside the powder bed while thermal energy provided through the newly deposited layer increases the part temperature in accordance with the heat equation given as where T is temperature, t is time, Q v is a volumetric heat source and ρ, c p , k are temperature-dependent density, specific heat and thermal conductivity, respectively. Next, to formulate a boundary value problem (BVP), process-relevant boundary and initial conditions are specified. It is reported that the effective conductivity of powder layer is only 1% of bulk conductivity [30,31]. Hence, it is justified to assume that the heat transfer through powder is negligible as compared to conduction within the part. Therefore, lateral sides of the part, i.e., the part-powder interface, is assumed to be thermally insulated where n i are the components of the outward normal unit vector to ∂Ω lat and the repeated index indicates summation over i. The part remains bonded to the baseplate which is typically pre-heated and acts as a heat sink. This is taken into account by specifying a temperature boundary condition where T 0 is the pre-heat temperature which is assumed to be constant. Next, heat losses through convection and radiation are considered at the top surface ∂Ω top where h conv is the convective heat transfer coefficient and T a is the ambient temperature, assumed to be constant within the machine chamber. The radiative heat transfer boundary condition is given as where σ is Stefan-Boltzmann constant and is emissivity of the radiating surface. Lastly, each newly added layer is assumed to be at an initial temperature of T 0 .
x 1 x 3 x 2 Figure 1. Schematic illustration of a body Ω being fabricated during the LPBF process. The topmost crimson colored region signifies the newly deposited layer. The body is attached to the baseplate at the bottom surface ∂Ω bot while a laser scans the top surface ∂Ω top . The part-powder interface is denoted as the lateral surface ∂Ω lat . Thermal losses due to convection and radiation are shown as q conv and q rad , respectively.
It is important to note that during the process, melting and re-solidification of material also takes place. However, the influence of these phase transformations on resulting temperatures is not considered here. This is motivated by the fact that the energy absorbed during melting is released during re-solidification when the laser beam moves away from the melt zone. Hence, the net effect of phase transformations is negligible when considering part-scale temperatures [32]. As discussed in Section 1, lumped layer-by-layer material addition is considered here, which means that the uniform heat source Q v is simultaneously applied over the entire volume of the newly deposited lumped layer. Both of these simplifications, i.e., lumping and layer-by-layer material addition, have been motivated and applied in number of previous LPBF simulation studies [23][24][25][26]33]. Apart from obvious computational benefits, these choices are motivated by the fact that it has been shown by Davies [24] that a model with these simplifications is adequate for correctly capturing the part-scale thermal response.

Finite Element Analysis and Simulation Parameter Setting
To solve the boundary value problem (BVP) given by Equations (1)-(5), finite element analysis is used. The element birth-and-death method is implemented which is a common technique used to simulate the growing domain during the AM processes [20,26]. The implementation is done in CalculiX (a free and open source FE analysis code) [34] with a structured mesh where elements are aligned with the (lumped) layers. Eight-node linear cubic elements are used where mesh size is equal to the superlayer thickness which is assumed to be S = 500 µm. This choice is motivated by the study conducted by Harrison [26] where temperature histories for different degree of layer lumping are compared. It is shown that the temperature response remains reasonably accurate for a superlayer thickness of 480 µm for Ti-6Al-4V parts with typical LPBF machine parameters. Typically, layer thickness for the LPBF process varies from 20-80 µm and in this study we use a representative value of 50 µm for real LPBF layer thickness i.e., the superlayer is 10 times thicker than the real layer thickness. The newly deposited superlayer is subjected to a uniform volumetric heat source Q v [W/m 3 ] for the heating time t h . In accordance with Equation (3), the finite element nodes associated with the bottom surface act as heat sink, i.e., their temperature remains fixed at T 0 . The heating step is followed by a cooling step in which the newly deposited layer can cool for an inter-layer time (ILT) t d . The cycle of sequential heating and cooling repeats while the thermal history of all previous heating/cooling steps influence the simulation of the next step.
As mentioned earlier, the modeling principles presented by Davies [24] and Davies [28] are used here which provides a basis for setting simulation parameters t h and Q v . The two key ideas presented by Davies [24] and Davies [28] are as follows: • Heating time t h for a layer with area A equals the time that the laser would take for scanning that entire layer, i.e., where h is hatch spacing and v is laser scan velocity. Please note that by using this definition, t h is calculated for each super layer and it depends on the local part geometry through layer scan areas.

•
It is ensured that the total deposited energy matches that of the actual process. Deposited energy per unit time is given by E = γP, where γ is the absorption coefficient and P is laser power. Using this principle, the volumetric heat source term can be calculated as where l is LPBF layer thickness.
Please note that the real LPBF layer thickness l which is different from the simulated superlayer thickness S is used for calculating Q v in Equation (7). This ensures that the input energy is automatically scaled when using thicker lumped layers. It is demonstrated by Davies [24] and Davies [28] that the described simulation scheme is capable of predicting the real physical temperatures recorded during the process as captured using the thermocouples located as close as 2.5 mm away from the topmost layer.
Next, the ILT parameter used for cooling down between layer depositions needs to be specified. As mentioned, S = 500 µm is used which is 10 times thicker than the real layer thickness. Hence, as suggested by Harrison [26], the ILT should also be scaled in order to compensate for the effect of large heat capacity of the superlayers. Hoelzle [23] used a method where ILT for a thick layer is adjusted such that thermal decay rates match with those of real-sized layers. However, the thermal decay rate varies at different locations within a part geometry which makes this estimation design dependent. Moreover, the ILT for the real size layers also depends on multiple factors such as layer area, scanning pattern, number of parts printed together inside the same chamber and recoater speed. Therefore, as per the linear scaling suggested by Harrison [26] and Malmelöv et al. [25], we choose to use an estimated value of 100 s which is 10 times the typical recoater time of 10 s and later we discuss the implications of this choice. The temperature-dependent properties of Ti-6Al-4V are taken from Davies [24] covering the range from room to fusion temperatures and shown in Figure 2. Finally, Table 1 lists the modeling parameters which are common to all the shown results.   Figure 3 illustrates the considered part geometry where dimensions and build direction are specified. The part dimensions are representative of a typical LPBF part. This particular design is chosen because of two key aspects. First, all overhanging features have the same overhang angle of 45 • , marked with red lines in Figure 3b. In fact, the design is obtained by topology optimization method with overhang angle control as presented by Langelaar [35]. For our purpose, we choose this design to investigate whether the same overhang angle exhibits similar overheating behavior or not throughout the sample. Secondly, the design features do not change in the depth direction, i.e., the design is simply an extrusion of the 2D design shown in Figure 3b. This allows for convenient visualizations of temperatures fields as temperatures do not change in the direction of depth when considering layer-by-layer simulations.   new layer addition is followed by a heating and cooling step, where the peak temperatures occur at the end of the heating step. As explained earlier, heating time t h varies as per Equation (6), while ILT t d remains constant. Please note that in this model the thermal history of all previous heating/cooling steps is stored and influences the next step. Figure 5 shows typical variation of the temperature with respect to time for nodes located at two different locations marked as A and B in Figure 3b. First, note that the temperature first rises and reaches to a maximum when a layer is activated and then it rapidly drops during the inter-layer time. Next, the temperature rises again when the next layer is added. This phenomenon repeats until the final layer is deposited. Secondly, note that peak temperature attained at Point B is significantly higher than that at Point A. This is due to the geometric features in the vicinity of Point B, which do not allow for efficient heat evacuation, hence resulting in a higher peak temperature. Also, Point A is located closer to the baseplate facilitating quicker heat evacuation. This suggests that peak temperature value can be used as an indicator of overheating risk associated with a geometrical feature. Figure 6 shows the peak temperatures attained for the entire domain during the build process. Please note that peak temperatures for different spatial locations may occur at different time instances during the simulation. It is evident that higher peak temperatures are exhibited near the overhangs, indicating local heat accumulation. Moreover, the funnel-shaped geometries in the region labelled as D in Figure 6 exhibit higher maximum temperatures than region labelled as C. This is in line with the experimental observations reported in Zimm. [17] and Patel et al. [18] where LBPF specimens exhibit local overheating in the vicinity of similar funnel-shaped features which act as a thermal bottleneck. This map of peak temperatures is used for judging design features for their heat accumulation behavior and referred to as the hotspot map. Moreover, although all features of the part have the same overhang angle of 45 • , the thermal behavior in their vicinity is not similar. This implies that use of a purely geometric design rule can be insufficient for avoiding local overheating.  It is important to note that in the actual AM process, the temperatures rise till the melting point and then stabilize due to phase change with excess heat resulting in a larger melt-pool. Therefore, some LPBF simulation studies put an upper limit on the temperature to address this phenomenon, e.g., Childs [21]. However, for the reference model we choose not to use this concept and instead use the overshoot of temperature as an indicator of heat accumulation tendency of the neighboring geometry. It is still important to mention that the temperatures found using this model are representative values and are not exact in-process temperatures, due to the simplifications already introduced in the presented reference model. However, the relevance of the presented model for judging design manufacturability is demonstrated by the presented example where it is shown that peak temperatures detect overheating tendencies associated with geometric features. Finally, the computational cost associated with the presented reference model remains significantly high due to the ever-growing analysis domain and the time integration to account for the transient nature of the problem. The model presented here is discretized using 2.16 million elements leading to 2.22 million nodal degrees of freedom (DOF) and the corresponding simulation time is 21 h 28 min 32 s. All computations reported in this paper are done on a HPC cluster node with 20 cores. The computational burden is still prohibitive for high-volume parts, interactive design iterations or integration with design optimization techniques. Therefore, to further reduce computational cost, additional physics-based simplifications are proposed in the next section.

Thermal Modeling Simplifications and Comparison Metrics
In total, six physics-based simplifications are presented in this section which are applied in addition to the simplifications already existing in the reference model detailed in Section 2. Section 3.1 discusses the influence of neglecting convective and radiative heat losses. Also, the effect of neglecting temperature-dependent properties is investigated. Although these simplifications have been applied in the literature, their impact in the context of detecting overheating was never studied. In Section 3.2, an analytical solution for one-dimensional heat equation is presented which serves as a basis for introducing three novel simplifications. With each simplification, a slightly different hotspot temperature field T sim is obtained. Hence, to compare these with the reference hotspot temperature field T ref , three comparison metrics are defined in Section 3.3.

Influence of Neglecting Convective and Radiative Heat Losses
As discussed in Section 2.1 and described by Equations (4) and (5), heat is lost through convection and radiation from the topmost layer. However, the relative importance of heat losses through convection and radiation as compared to the conduction within the part is still debated in the literature. Some studies present arguments in favor of neglecting these losses [19,20,23], while others advocate for their inclusion in the thermal modeling of the AM processes [24,36]. A basic estimate can be made to quantify the thermal losses due to each mode of heat transfer. The convective loss near the melt zone can be estimated as where T m is the melting point of Ti-6Al-4V taken as 1604 • C and other parameters are given in Table 1. Similarly, radiative loss can be estimated as Next, in order to compare, the rate of energy transfer within the part due to conduction, q cond can be estimated using a simplified version of Fourier's law for one-dimensional heat flow: where ∆T is the temperature difference measured over a distance ∆x. The melt zone is considered for making this estimation as the highest thermal gradients occur there. Consequently, conductivity at the melting point is considered in Equation (10). Due to extremely high thermal gradients and optical inaccessibility, it is extremely challenging to accurately record actual temperatures near the melt zone. Nevertheless, available data from the literature can be used to make an estimate. For example, empirical observations from Davies [28] reported a temperature difference ∆T = 1200 • C between the topmost point of the melt zone and a point located 2.5 mm below it, i.e., ∆x = 2.5 mm. Putting these values in Equation (10) gives |q cond | = 1.4 × 10 7 W/m 2 . Please note that these are only rough estimates as heat fluxes continuously change with time-varying temperature fields. Moreover, this estimation of q cond is an underestimation as even higher temperature gradients are observed in the vicinity of melt zone by melt-pool simulation studies [37]. Nevertheless, it can be observed from these estimates that indeed heat losses through convection and radiation are orders of magnitude lower near the melt zone than those by conduction, as argued by Paul et al. [19]. Also, radiation accounts for more heat loss than convection. The effect of neglecting one or both of these loss terms will be studied.
In the remainder of this paper, simplification by exclusion of radiation is referred to as S1 and that for excluding convection is termed as S2.
Another aspect that makes the PDE given by Equation (1) nonlinear, and consequently computationally expensive, is the temperature dependence of the thermal properties. Hence, another simplification (S3) is considered where along with exclusion of convection and radiation, constant thermal properties are considered. Ayas [32] and Yang et al. [38] performed a calibration study and showed that use of constant melting point properties are suitable when probing local temperatures near the heat deposition zone. Hence, in simplification S3, constant values of ρ = 4200 kg/m 3 , c p = 750 J/kg K and k = 27.5 W/m K are used. Results for detecting hotspots and computational gains achieved using these simplifications are reported in Section 4.

Novel Simplifications Motivated by One-Dimensional Heat Transfer Analysis
To gain fundamental insights in the nature of transient heat transfer phenomena relevant to the LPBF process, a simplified case of one-dimensional heat transfer is considered and the analytical solution is presented for the boundary conditions analogous to the reference model. The detailed derivation is presented in Appendix A while the final solution and its implications are discussed here. A one-dimensional domain with length L is shown in Figure 7a. The heat equation given by Equation (1) can be simplified for one-dimensional heat transfer as where x is position and α is the thermal diffusivity given by α = k/ρc p . Please note that the temperature dependence of thermal properties has been neglected for simplicity and constant values as given in Section 3.1 are used. First, a heating step is considered in which the topmost point x = L is subjected to a heat flux Q analogous to the volumetric heat source used in the reference model. Next, a heat sink condition is assumed at the bottom, i.e., x = 0. To derive a more general solution, the sink temperature is specified as T s unlike the reference model where it was assumed same as the initial temperature T 0 . These boundary conditions are given as The initial condition is x L Q

Time
Temperature The rod is heated for a time t h followed by a cooling step in which Q = 0. Using the method of separation of variables, the analytical solution is derived in Appendix A reads where T h represents the temperature distribution during heating step and T h (L, ∞) = QL/k represents the steady-state temperature value at x = L while λ n = (2n − 1)π/2. Next, the rod is allowed to cool, i.e., the boundary condition given by Equation (12) becomes while the other boundary condition remains the same as given by Equation (13). For this case, the temperature distribution at the end of the heating step becomes the initial condition, i.e., Again, the PDE given by Equation (11) is solved and the solution for the cooling regime is given as where T c represents the temperature distribution during the cooling step. A visual depiction of the derived equations is given in Figure 7b where temperature variation for the topmost point (x = L) with respect to time is shown. Below, the derived solutions will be used to motivate three novel simplifications, aimed at further reduction of the computational burden associated with the reference model.

Observation 1: Temporal Decoupling
First, we focus our attention on the magnitude of the temperature drop that occurs in the cooling phase between layer depositions, using the one-dimensional model. When this drop is sufficiently large, a decoupling in time of deposition steps can be considered in the process model. Recall from Section 2 that the peak temperatures at the end of the heating step are used to construct the hotspot map. Hence, the temperature difference between x = L and x = 0 is considered. This difference during the cooling regime is normalized with the maximum temperature difference attained at the end of the heating step. The normalized temperature differenceT c (L, t) for the topmost point then readŝ Please note that the sink temperature is assumed to be the same as the initial temperature in accordance to the boundary conditions used for defining the reference model, i.e., T s = T 0 . It can be deduced from Equation (19) that the cooling behavior depends on the duration of the heating stage t h and the characteristic time of the heat equation τ = L 2 /α. Figure 8a,b show the variation ofT c for a range of t h and τ values, respectively. These ranges are selected in the context of the LPBF process. For example, t h is typically very short in LPBF considering a typical scanning velocity v = 1 m/s with a laser spot diameter of 100 µm. Therefore, t h = 10 −3 , 10 −2 , 10 −1 s are selected to demonstrate the effect of heating time. Similarly, parts as high as 300 mm can be built in LPBF machines giving τ max ≈ 10 4 s, assuming α for Ti-6Al-4V. Hence, this value is used along with lower values to study the variation. Also, constant τ = 10 4 s and t h = 0.01 s are used for plotting graphs for varying t h and τ, i.e., Figure 8a,b, respectively. The infinite series in Equation (19) is converging and n = 10 4 is found to be sufficient. Both graphs show that a slower cooling is observed as t h or τ increases. The graphs are shown for the first 10 s of cooling which is close to typical ILT values used in LPBF [11,23,26]. The crucial observation here is that in both graphs the topmost point cools down to approximately 10-20% of its highest temperature value in this time frame. This suggests that topmost point cools down to an estimated value of 150-300 • C assuming T m as the highest temperature.
It is noteworthy that this finding is based on the simple one-dimensional model which assumes a constant conductivity throughout the domain. However, in the real three-dimensional setting, there could be local zones of lower conductivity near the topmost layer which would add to the heat accumulation and decrease the local cooling rate. As described in Section 1, examples of such features are acute overhangs and thin connections. Nevertheless, it is found that the findings based on this one-dimensional model are in agreement with the experimental observations reported by Davies [28] where thermocouples embedded inside an overhanging part recorded temperatures during the build. It is shown that the part temperatures remain in the range of 100-400 • C at different locations below the topmost layer when the pre-heated baseplate temperature is at 100 • C. This suggests that the one-dimensional insights can be extended to real parts. Another interesting observation made on the same data is that the recorded temperatures increase with the height at which thermocouples are placed. The effect of build height is also investigated by Hilgenberg [11] and slower cooling is reported for increasing part height. As a qualitative comparison, this is in accordance with cooling curves shown in Figure 8b, since τ ∝ L 2 . It suggests that the effect of build height and thermal properties can be combined as characteristic time τ in order to study the cooling behavior. A more detailed quantitative comparison between τ and cooling data will better establish this correlation. However, it is deemed out of the scope of this paper. τ = 1 s τ = 10 s τ = 10 3 s τ = 10 4 s The observations from the one-dimensional analytical model and experimental data from literature indicates that the ILT in the LPBF process allows for ample cooling between the layers. Moreover, it is also recommended as good manufacturing practice to design process such that appropriate cooling happens between two successive layers [11]. This suggests that the previously deposited layers do not significantly affect the peak temperatures recorded for the next layer. This enables decoupling the thermal history of different layers from each other for peak temperature prediction, which essentially means that each layer addition can be assumed to have an initial temperature equal to the baseplate temperature T 0 . Peak temperatures attained this way would still capture the local overheating associated with the geometrical features of the part. A schematic representation of this idea is given in Figure 9a, where previous thermal history is not transferred to the next layer addition. There are two major computational benefits associated with this simplification. First, only simulation of the heating step suffices for capturing the first peak temperatures for each layer. Second, heating steps associated with every layer of the structure can be computed in parallel. The model which makes use of this simplification is referred to as the 'temporally decoupled model' and represented as S4. It is important to note that this model cannot capture the gradual heat accumulation that may occur over layer depositions, as information about the thermal history is lost. Nevertheless, it is found that features prone to local overheating can still be quickly and adequately identified when making use of this simplification. Recall that these features also contribute significantly to the gradual heat build-up that happens over the layers.

Observation 2: Spatial Decoupling
Another useful observation is made by focusing on the relationship between the peak temperature at the end of the heating step T h (L, t h ) and domain size L. This implies substituting t = t h , x = L and T s = T 0 in Equation (15) which gives This relationship is pictorially presented in Figure 10 for three different heating times. It is evident that there exists a saturation domain length L s and increasing the domain size beyond this value has no effect on the peak temperatures. This phenomenon is characterized by the exponent term in Equation (20) which constitutes the Fourier number given as Fo = αt h /L 2 . Please note that Fo is reducing along the horizontal axis in Figure 10 as L increases. It can be deduced from the graph that the saturation regime starts at Fo = 0.3 which is marked by the vertical lines. Physically, it implies that for a given domain size L, thermal diffusivity α and heating time t h , if Fo ≤ 0.3 then considering larger spatial domains will not influence the peak temperatures. It follows that the corresponding saturation length is L s = 1.82 √ αt h .

L [mm]
T h (L, Figure 10. Variation of peak temperatures attained at the end of the heating step with respect to domain size L as described by Equation (20). Three graphs are shown for varying heating time t h and vertical lines are shown for respective Fo = 0.3. First 10, 000 terms of the infinite series given by Equation (20) are considered for plotting.
In the context of LPBF modeling, this observation implies that during the heating step, if the domain beyond L s is discarded from the analysis, peak temperatures will not be affected. Although derived using the one-dimensional model, the idea can be easily extended to higher dimensions where low conductivity regions influence the peak temperatures, only if they are present in the vicinity of topmost layer. Recall that temperature dependence is neglected, and the value of thermal properties estimated at the melting point are used for analytical derivation. However, since the highest value of α is achieved at the melting point, this choice gives an upper limit for the saturation length. Figure 9b shows an implementation of this idea, where the Dirichlet boundary condition given by Equation (3) is now applied at a distance L s below the heat source instead of at the baseplate. This simplification will be applied together with the previously introduced simplification of decoupling the layers in time. Hence, it is referred to as 'spatially and temporally decoupled model' and represented as S5. The reduced domain size enables further reduction of the computational cost.

Observation 3: Steady-State Response for Detecting Overheating
It is well-known from the theory of heat conduction that lower conductivity regions inside a domain would influence the steady-state response of a thermal analysis [39]. Consequently, a steady-state response can also be used for detecting regions of lower conductivity within a given domain. The computation of steady-state responses is much faster than a transient analysis which makes it an attractive option for quickly detecting regions prone to overheating. In the context of the one-dimensional rod, substituting t = ∞ in Equation (15), the steady-state temperature distribution along the length of the rod is found as This linear thermal profile is depicted in Figure 11a where x varies from 0 to L along the vertical axis while temperatures are plotted along the horizontal axis. Another case shown in Figure 11b considers that a subsection of the one-dimensional rod has lower thermal conductivity than the rest of the domain. This is shown by the orange patch at Location A. The steady-state temperature distribution for this case is computed and plotted, where a higher thermal gradient is observed in the patch of lower conductivity. Please note that also a higher steady-state temperature at x = L is found compared to the case without the patch, indicated by the dotted line. This signifies that the steady-state temperatures can also be used for identifying regions with lower conductivity which are prone to overheating. In three-dimensional setting, examples of such features are acute overhangs and thin connections which would create a local zone of lower conductivity. However, the steady-state response should be cautiously used as it does not take into account the proximity of low conductivity regions to the topmost point where the heat flux is applied. This is illustrated in Figure 11c where the same low conductivity patch is now situated at Location B and still causes the same increase in the steady-state temperature at x = L. This implies that a low conductivity region situated far away from the heat deposition zone has the same effect on the top temperature when situated close to the topmost point, when analyzed using steady-state response. In context of the three-dimensional parts, this is unrealistic as heat transfer phenomenon is transient and only nearby regions of low conductivity (e.g., acute overhangs, thin connections) would influence the thermal behavior near the melt zone. Therefore, in order to rectify this inherent limitation of steady-state analysis, only the domain close to the topmost layer must be considered. Hence, this simplification is applied together with previously explained simplifications of temporal and spatial decoupling. In other words, instead of performing a transient analysis in the spatially and temporally decoupled model shown in Figure 9b, a computationally fast steady-state analysis is performed and the simplified model which involves this approximation is referred to as 'steady-state model' represented as S6.

Comparison Metrics
To assess the ability of each proposed model to detect hotspots in a part, three different comparison metrics are defined in order to capture different aspects associated with overheating detection. Recall that reference and simplified temperature fields are referred to as T ref and T sim , respectively. The first metric is defined as the percentage error between the maximum temperatures recorded for both fields, This quantity compares the intensity of worst overheating as predicted by the reference and the simplified models. It indicates how much the corresponding simplification under/over-predicts as compared to the reference case. A positive δ signifies over-prediction and vice versa. Please note that this quantity compares the absolute maxima between two temperature fields without requiring that these maxima occur at the same spatial point. Therefore, to assess the spatial similarity between two temperature fields, the Jaccard index is used. This is a measure commonly used in the field of image recognition for quantifying similarity between two images and is defined as the ratio of intersection and union between the two images [40]. For our case, first T is normalized using max(T), i.e.,T = T/max(T), whereT indicates the normalized temperature field. Next, the Jaccard index is defined as where T sim represent normalized temperature values for node i for the reference and the simplified fields, respectively and n is the total number of nodes in the finite element analysis. A high Jaccard index implies that the overall temperature distributions over the part are highly similar.
The third metric is defined for a qualitative comparison of two hotspot fields. A critical zone identification (CZI) map is defined which highlights the worst overheated zones while suppressing the cooler, less critical, regions in a given geometry. This allows for judging simplifications based on their capability of detecting the correct worst overheating locations while the actual temperature predictions might be significantly off. It is defined as the contour map of the normalized temperature field with four contour levels, here chosen at 50%, 70%, 80% and 90% of the maximum temperature in the field. Using this, different design features can be identified from least to most critical in terms of overheating.
Finally, in order to quantify the degree of computational gain achieved by a simplification, a gain factor is defined as where t c ref and t c sim are the wall-clock times for completing the reference and simplified analysis, respectively.

Numerical Results and Discussion
To illustrate the validity of the proposed simplifications and discuss the associated implications, hotspot maps for the part shown in Figure 3 are generated using simplified models. For ease of visual comparison, all hotspot fields are presented in Figure 12 along with corresponding CZI maps. Simulation times and calculated comparison metrics are listed in Table 2 while the default simulation parameters are used from Table 1. All computations are performed using 20 cores on a HPC cluster. Discussion about each simplification is presented in the same order in which simplifications were introduced in Section 3.

Hotspot Map without Considering Convective/Radiative Heat Losses
To compare the relative importance of convection and radiation heat losses, two hotspot maps are prepared. One by excluding radiation, shown in Figure 12c and other by excluding convection, shown in Figure 12e. All other simulation parameters remain the same as in the case of the reference model, depicted in Figure 12a. The first observation is that neglecting convection or radiation leads to a conservative estimation of overheating zones. This is due to the fact that both convection and radiation contribute to the thermal losses and hence lead to higher simulation temperatures, when excluded. Also, as estimated, exclusion of radiation has a greater impact than exclusion of convection. It is manifested by the fact that the maximum temperature attained by excluding radiation is 22.4% higher than that obtained with the reference model, i.e., δ(S1) = 22.4%. On the other hand, δ(S2) = 4.2% when convection is excluded from the analysis. Moreover, when looking at the overall temperature distribution over the entire domain, the hotspot map obtained considering radiation while neglecting convection is very similar to that obtained using reference model. This is supported by the Jaccard index J(S2) = 96.4%. Please note that even with exclusion of convection or radiation, the hotspot map still correctly predicts the funnel-shaped geometries as the worst overheating features. This is also clear from the CZI maps shown in Figure 12b,d,f where the same features are identified as critical. Lastly, excluding radiation and convection provides marginal computational gains of η(S1) = 1.06 and η(S1) = 1.03, respectively. The radiation boundary condition is nonlinear and hence its exclusion provides a slightly higher computational gain.

Hotspot Map Without Considering Temperature-Dependent Properties
The influence of neglecting temperature dependence of thermal properties is studied next. The hotspot map presented in Figure 12g considers this simplification in addition to the simplifications of excluding convection and radiation. The thermal properties as listed in Section 3.1 are used for preparing this map. An interesting observation here is that a significant heat accumulation is observed between the layers. This is indicated by the temperature gradient that is seen in the build direction where the topmost layers exhibit significantly higher peak temperatures than the reference model. Due to this, the overall distribution of the peak temperatures over the domain is considerably different from those found using the reference model. This difference is quantified by J(S3) = 75.4%. Moreover, the maximum temperature is overestimated by δ(S3) = 19%. Nevertheless, the funnel-shaped geometries are still detected as the worst overheating zones, as shown in the CZI map presented in Figure 12h. The computational gain factor is found to be η(S3) = 2.1.
The observations made using these simplifications can be used for making appropriate modeling choices when investigating different aspects of the LPBF process. It is shown here that the exclusion of convection, radiation and temperature dependence halves the simulation time while correctly detecting design features responsible for severe local overheating. However, if the aim of the simulation is to study the effect of gradual heat accumulation during the build, for example as conducted by Jamshidinia and Kovacevic [41], then using these simplifications might lead to overly conservative fields and, thus, leading to false positives. This will become more critical in the case of tall parts with high heat capacity and/or processes with short ILT. Therefore, suitable modeling choices should be made.

Hotspot Map with Temporal Decoupling
Inspired by the observations from the analytical model, the simplification of decoupling layer addition is studied. The hotspot map in Figure 12i shows the peak temperatures obtained at the end of the heating step. The computations associated with each intermediate layer addition are carried out in parallel and the hotspot map is prepared as a post-simulation step. As discussed in Section 3, this simplification assumes constant initial temperature T 0 for each new layer addition. In other words, it assumes that the part cools down to T 0 between successive layer additions. Consequently, this simplification cannot capture the effect of gradual heat accumulation which builds up over the layers. In fact, this effect is opposite to the preceding case (S3) where gradual heat build-up was overestimated. This is also evident by the fact that peak temperatures are underestimated here, as quantified by δ(S4) = −5.8%. Nevertheless, the hotspot map prepared using this simplification yields result very similar to the reference case. This is evident by comparing the CZI maps shown in Figure 12b,j which highlight same set of design features. This is also quantified by high similarity between the two temperature fields with J(S4) = 89.9%. Lastly, a computational gain factor η(S4) = 85.2 is achieved where thermal analysis takes only 15 m 7 s while the difference between peak temperature values δ(S4) remains less than 6%. Recall that reduction in wall-clock time is attributed to the parallel simulation of each new layer addition and omission of the cooling step simulation between the layers. As mentioned, this model cannot capture the gradual heat accumulated during the build which is inversely proportional to the ILT [11], and similar precautions apply as mentioned for S3.

Hotspot Map with Spatial Decoupling
As per the observation discussed in Section 3.2.2, the simplification of considering only the thermally relevant domain during the heat addition step is applied and the resulting hotspot map is presented in Figure 12k. Recall that this simplification is applied in addition to the simplification of decoupling in time, as shown in Figure 9b. The domain size for each new layer addition is calculated using the concept of saturation length presented in Section 3.2.2. The resulting hotspot map is found to be very similar to the one found using the concept of temporal decoupling (S4). This signifies that very small additional errors are introduced when considering local domains instead of the full geometry. Reduction of the size of simulation domain provides a further improvement in computational gain factor, which reaches 144.2.

Hotspot Map with Steady-State Analysis
Next, use of computationally fast steady-state analysis instead of the transient response is studied. The hotspot map shown in Figure 12m shows the obtained hotspot map with steady-state peak temperatures. As expected, the temperature values are significantly higher than those found using the transient analysis. Moreover, the overall temperature distribution is also considerably different. These are quantified by δ(S6) = 65.3% and J(S6) = 74.8%. Nevertheless, the hotspot map found using steady-state analysis also identifies the correct critical zones of overheating, as observed in the CZI map presented in Figure 12n. The analysis is completed in only 2 min 9 s with a computational gain factor of 600.

Comparative Analysis
A pictorial representation of maximum percentage error δ and Jaccard index J obtained using different simplifications is given in Figure 13a and computational gain factors are plotted in Figure 13b. Please note that δ is a measure of difference between two fields while J is a measure of similarity. Therefore, the simplifications which yield low δ and high J signify higher conformance with the reference hotspot map. Upon comparison, it becomes clear from Figure 13a that the analysis which includes radiation (S2) is more accurate than the one which excludes it (S1). Also, the analysis with decoupled layers (S4) and local domain (S5) are almost the same in terms of accuracy. Figure 13b highlights the considerable computational gains provided by the novel simplifications proposed in this paper. Please note that these gains are achieved in part due to the parallel processing, which becomes possible due to the proposed simplifications. The total processing CPU times are also reported in Table 2 which presents the total computation time used by all the processors. This also implies that the wall-clock time will directly depend on the number of processors used. Please note that even without using parallel processing, the novel simplifications provide significant gains where the steady-state model is still 15 times faster than the reference model when CPU times are compared. To conclude, a summary of disadvantages or risk associated with each simplification is presented in Table 3.
Percentage value Computational gain factor |δ| J

Conclusions and Future Work
In this paper, first an established thermal model is used to predict the overheating behavior in an LPBF part. A representative part with a 162 cm 3 deposition volume is considered which is simulated using 2.22 million DOFs. It is revealed that using purely geometric design guidelines might not be sufficient for avoiding overheating. It is shown that even with the simplifications of layer lumping and simultaneous simulation of entire layer deposition, the computational time is still prohibitive for quick manufacturability assessment needed for efficient design modifications, process-parameter tuning and optimizations. Hence, a total of six further simplifications are investigated, ranging from omission of radiation and convection and use of constant material properties, to novel simplifications involving temporal and spatial decoupling and ultimately the use of localized steady-state responses. It is shown that in particular the three novel simplifications provide very high computational gains. The results from the simplifications of temporal and spatial decoupling show that even after omitting the simulation of the cooling step and reducing the computational domain size, these simplifications can still provide crucial information about design features and their thermal behavior during the LPBF process. Using these simplifications, the correct locations prone to overheating can be identified and error in maximum temperature prediction is less than 10% when compared with the reference case. Moreover, the localized steady-state response is shown to provide accurate qualitative information about the location of problematic features while it is 600 times faster than the reference model. Note that this computational gain directly depends on the number of processors used.
The high computational gains achieved using simplified models makes them especially suitable for design optimization problems where several hundred design evaluations might be needed for finding the optimized design. For example, integration of simplified models with topology optimization techniques holds the promise to deliver highly efficient designs which are also robust from the context of overheating. Typically, design optimization methods make use of gradient information and it has been shown by Van Keulen et al. [42] that gradient computation for a transient analysis is computationally much more expensive as compared to a steady-state analysis. Consequently, the steady-state model which correctly identifies overheating zones is a perfect candidate for integration with topology optimization, which is seen as the immediate next step.
The layer-by-layer model used in this research cannot capture the influence of laser scanning vectors on local overheating. Therefore, development of higher fidelity models which simulate laser movement and analysis of the generated thermal history for identification of local overheating is seen as a future step. It is expected that the introduced simplifications of temporal and spatial decoupling would remain valid for such higher fidelity models. However, thorough investigation of their validity would still be needed and seen as a possible research step. The thermal history also dictates the development of residual stresses which cause deformations and even part failures. Therefore, the possibility of extending simplified thermal models for capturing mechanical aspects is a promising research direction. Another avenue of future research is to experimentally validate the simplified models. For this purpose, a study is underway where the beam design investigated in this paper is manufactured and a metallographic study is being performed. The specimen will be cut and examined near the funnel shapes which are identified as critical zones by the simplified models. Now we apply principle of superposition, which means that if the solution holds true for one λ, then it must be true for a linear combination of all λ n . Hence, Finally, we substitute D n from Equation (A19) into Equation (A13) and subsequently use Equation (A3), to get temperature distribution during the heating regime as Note that T ss = (QL)/k is same as the steady-state temperature at x = L for the heating step, i.e., T ss = T h (L, ∞).
The derivation for cooling regime remains exactly the same except the initial condition is now prescribed by the temperature distribution at the end of the heating step, i.e., T(x, 0) = p h (x) = T h (x,t h ) forx > 0.
(A21) Therefore, replacing p(x) by p h (x) in Equation (A17), integrating and using Equation (A3) leads to the solution for cooling regime as