Introduction

Laser powder bed fusion (L-PBF), which is one of the metal additive manufacturing (AM) processes, can fabricate near-net shape products with high precision. In the L-PBF process, a thin layer of raw metal powder on a platform is scanned by a focused high-power laser beam to selectively melt and solidify the powder. Powder spreading and laser beam scanning are repeated until the fabrication is completed. The advantages of L-PBF in the aerospace and aircraft industries include significant cost and lead-time reductions, novel material development, unique design solutions for efficiency and lightweight, and consolidation of multiple components for performance enhancement and risk management [1]. For such practical applications, it is essential to understand complex relationships between process, structure, and property (P–S–P) resulting from the interaction of multiple physical phenomena at different spatial and temporal scales during the process. The P–S–P relationships for some typical metals (e.g., alloys of titanium, nickel, and aluminum) have been well studied using experimental approaches and simulations over the past decade, so that the next targets are to improve product performance by process parameter optimization [2,3,4,5,6,7,8,9,10] and exploration of novel materials for the process [11]. Since such optimization and exploration strategies require many trial-and-error iterations, experimental approaches are expensive and time-consuming. Thus, it is necessary to use physics-based modeling on computers that reproduces the L-PBF process [2, 9, 12,13,14] coupling with subsequent simulations, such as thermo-elastoplastic analysis for solidification cracking [9], thermo-mechanical analysis for residual stress distribution [2], and microstructure evolution [12].

The L-PBF process is usually modeled by coupling computational fluid dynamics (CFD) with powder spreading simulations and ray tracing techniques for laser beam absorption [12,13,14]. Such thermo-fluid analysis reveals the accurate melting pool geometry with temperature distribution as well as lack-of-fusion [12], keyhole formation [13], and spattering [14]. On the other hand, such CFD-based simulation generally takes a considerable amount of computational time even with parallelization. In contrast, finite element method (FEM)-based thermal analysis with a moving heat source model is more computationally efficient and suitable for iterative procedures such as process optimization [9, 15].

Instead of strictly modeling the fluidity of melting pool and laser beam absorption based on actual phenomena in CFD-based simulation, a heat source model for the FEM-based thermal analysis is simply defined by a volumetric heat flux with effective absorptivity and shape parameters. For years, various heat source models with different heat flux distributions have been proposed and used for metal AM processes as well as welding processes. In 1946, Rosenthal was the first to derive an analytical solution for a moving heat source to roughly understand the temperature field in solids [16]. The solution was applied to obtain the approximate width of the melting pool in the L-PBF process [17, 18]. However, the analytical solutions do not provide the correct temperature field with consideration of temperature-dependent physical properties, boundary conditions, and heat source scanning strategies. Promoppatum et al. [18] have comprehensively compared the Rosenthal equation and the FEM-based thermal analysis with a two-dimensional Gaussian heat source model, and concluded that both approaches provided similar and reasonable thermal and microstructural results at a low-energy input, but the Rosenthal equation should be used with caution at a high-energy input. Unlike these analytical solutions, a double ellipsoidal heat source model for finite element thermal analyses was originally proposed by Goldak et al. [19] to represent both shallow and deep penetration of a melting pool in welding processes [19, 20], and is still most commonly used in the L-PBF process [21,22,23]. Even though there are several other heat source models (such as models with conical [24, 25] and depth-decaying geometries [7, 26,27,28], a hybrid model for the conduction and keyhole modes [29]), it is not fully understood which model is the most appropriate for the L-PBF process. This is because an appropriate calibration method has not been well established, and heat source models have not been fairly compared based on a quantitative indicator. Zhang et al. [30] performed finite element thermal analysis with eight different heat source models, but all eight models resulted in melt pools that were 40% shallower than in the experiment. This may have been due to the heat source models not being properly calibrated.

In order to validate the thermal analysis, effective absorptivity and shape parameters in a heat source model should be calibrated by comparing the simulation and experimental results. The effective absorptivity in the L-PBF process is higher than the absorptivity at a flat bulk metal surface because the laser beam is reflected multiple times in the powder layer and at the surface of a concave melting pool caused by recoil pressure [14, 31]. One method to directly measure the effective absorptivity is to use calorimetric measurement with thermocouples under the powder layer, as established by Rubenchik et al. [31,32,33]. The measured absorptivity of the powder layer was about twice as large as that of the flat surface. Using the same measurement method, Trapp et al. [31] found that the absorptivity was significantly dependent on laser power and scanning velocity (e.g., from 0.35 to 0.7 for 316L stainless steel) due to the transition from conduction to keyhole mode. While such measurements require a specially designed platform, the current study proposed a simplified approach based on the above analytical solution and an experimental bead-on-plate test.

On the other hand, the shape parameters were often calibrated by comparing the melt pool length, width, and depth by thermal analysis and a bead-on-plate test [20, 22,23,24,25,26,27, 29, 30, 34]. However, it has not been quantitatively determined whether the shape parameter values are the most appropriate. Ross et al. have introduced a systematic inverse heat conduction problem (IHCP) method to calibrate a volumetric heat source model to work for both conduction and transition melt pool modes [21]. Their study found clear relationships between absorptivity and shape parameters in the heat source model and the melt pool depth and length in the experiment. In our previous study [35], the simulation error was quantitatively defined by taking the difference of fusion zone geometry between thermal analysis and sample cross-sectional observation. Then, the optimal shape parameters to minimize the error were found by the iterative thermal analysis with Bayesian optimization. Due to this calibration strategy, the shape parameters with an error of less than 4% were successfully found with the minimum computational effort. The similar calibration strategy successfully found the optimal heat source model parameters for the welding process [20].

Such a calibration procedure should be performed for each laser scanning condition because the effective absorptivity and melt pool shape depend on the laser scanning conditions. As pointed out by Bayat et al. [36], the need for re-calibration limits the applicable window of laser scanning conditions in the finite element thermal analysis. This is especially crucial for optimizing the laser scanning conditions by the iterative thermal analysis because it is impractical to re-calibrate the heat source model parameters every time the conditions change. Thus, if the effective absorptivity and shape parameters are properly expressed as a function of the laser scanning conditions, the limitation will be eliminated and the thermal analysis will actually be valid with an arbitrary condition within the process window.

In this study, therefore, we propose a novel heat source model that can be applied to a wide process window of L-PBF. The secondary objective is to quantitatively discuss whether the heat source models developed to data are optimal for the L-PBF process. To these ends, we first conducted the bead-on-plate test using IN738LC as the model material (Sect. “Bead-on-plate test”). Then, the effective absorptivity was calibrated based on the Rosenthal analytical solution (Sect. “Estimation of effective absorptivity”). Also, the shape parameters were calibrated by the iterative thermal analysis with Bayesian optimization (Sect. “Screening heat source models”). These calibration procedures allowed to fairly and quantitatively determined the best of the four heat source models. Finally, the effective absorptivity and the shape parameters were represented as a function of laser power and scanning velocity using a multiple linear regression analysis (Sect. “Functionalization of shape parameters as laser scanning conditions”). Using the proposed heat source model, a case study was also demonstrated in Sect. “Discussion” to search for the optimal laser scanning condition for achieving both target fusion depth and cooling rate.

Materials and Methods

Bead-on-Plate Test

A bead-on-plate test for a nickel superalloy, Inconel 738LC (IN738LC) was conducted using a commercial L-PBF machine (SLM 280; SLM Solutions). A plate was placed on the platform, and powder was spread over it by the scraper. The powder layer thickness was roughly regarded as 30 µm, which is the thickness of one deposition layer in the L-PBF process. The IN738LC plate was prepared by cutting from an ingot cast in the laboratory. As for IN738LC powder, AMPERPRINT 0151.074 by Höganäs AB, with diameters ranging from 15 to 45 µm, was used. Table 1 shows the chemical composition of IN738LC provided by the manufacturer.

Table 1 Chemical composition of the IN738LC alloy (wt.%, Ni = balance)
Full size table

Single beads with a length of 10 mm were fabricated on the base plate by laser beam irradiation at different laser powers \(P\) and scanning velocities \(v\). As shown in Fig. 1b, the laser power \(P\) and scanning velocity \(v\) were varied from 100 to 700 W and 500–2000 mm/s, respectively, and the laser scanning conditions were selected such that \(P/v\) was between 0.133 and 0.5 J/mm. The wavelength and beam diameter of the laser were 1.064 and 80 µm, respectively. The chamber atmosphere was replaced with argon gas throughout the bead-on-plate test to keep the O2 level below 0.01 volume percent. The bead samples were cut in the middle of the length direction, and the polished cross section was observed by using scanning electron microscopes (SEM) (JSM-6010 LA and JSM-7200F; JEOL). As in the previous study, the cross-sectional SEM image was sequentially processed with the software Fiji [37] to crop the region of interest, binarize the fusion zone and the surrounding non-melt area, and make it coarse-grained to match the mesh size for the thermal analysis. As shown in Fig. 1a, these image processing created a two-dimensional array \({{\text{FZ}}}_{{\text{obs}}}\) with components 1 and 0 representing the fusion zone and the unmelt area, respectively. In addition, the area, width, and depth of the fusion zone were measured and assigned as \({A}_{{\text{obs}}}\), \({w}_{{\text{obs}}}\), and \({d}_{{\text{obs}}}\), respectively (see Fig. 1a).

Fig. 1
figure 1

a A cross-sectional SEM image of the bead-on-plate test was processed to prepare a two-dimensional array \({\text{FZ}}_{{{\text{obs}}}}\) with components 1 and 0 representing the fusion zone and the unmelt area, respectively. As indicated in the SEM image, area \(A_{{{\text{obs}}}}\), width \(w_{{{\text{obs}}}}\), and depth \(d_{{{\text{obs}}}}\) were also measured. b Laser scanning conditions for the bead-on-plate test. The color and size of the plotted points indicate the laser power and scanning velocity, respectively. The dark and light gray dotted lines represent \(P/v\) of 0.5 and 0.133 J/mm, respectively

Full size image

Analytical Solution

The following equation represents temperature T in a solid at time t after a point heat source moving passes through in the \(x\) direction, which was originally derived by Rosenthal [16]:

$$ \begin{array}{*{20}c} {T = T_{0} + \frac{\eta P}{{2\pi kr}}\exp \left( { - \frac{{v\left( {r + \xi } \right)}}{2\alpha }} \right)} \\ \end{array} $$
(1)

where \({T}_{0}\) is the temperature at locations far from the heat source, \(\eta \) is the effective absorptivity, \(k\) is the thermal conductivity, and \(\alpha \) is the thermal diffusivity. \(r\) is the distance from the heat source and defined as \(\sqrt {\xi^{2} + y^{2} + z^{2} }\), where \(\xi = x - vt\). It should be noted that this analytical solution was derived based on the following assumptions: (i) it is assumed to be a point heat source with constant \(P\) and\(v\), (ii) the material properties \(k\) and \(\alpha \) are temperature independent, and (iii) phase transformation, convection in a liquid (the melting pool), heat loss from surface convection, and radiation are neglected.

The analytical solution was used to calibrate the effective absorptivity\(\eta \). First, as shown in Fig. 2a, the analytical temperature field was simulated using Eq. (1). Then, the \(yz\) cross-sectional fusion zone where the temperature exceeded the liquidus temperature \({T}_{{\text{l}}}\) through the heat source passed was determined (see Fig. 2b). Qualitatively, it is obvious that the fusion zone area \({A}_{{\text{as}}}\) increases as \(\eta \) increases. Thus, the value of \(\eta \) at which \({A}_{{\text{as}}}\) is equal to \({A}_{{\text{obs}}}\) was determined by the bisection method. In other words, \(A_{{{\text{as}}}} \left( \eta \right)\) was repeatedly calculated until \(\left| {F\left( \eta \right)} \right| = \left| {A_{{{\text{as}}}} \left( \eta \right) - A_{{{\text{obs}}}} } \right|/A_{{{\text{obs}}}}\) was below \(\delta =\) 0.005. As in the previous studies for the melting pool width prediction with the analytical solution [17, 18], the material properties of \(k\) and \(\alpha \left( { = k/\rho C_{p} } \right)\) for IN738LC at room temperature \({T}_{0}\) were used for this calibration (see Sect. “Finite element thermal analysis” for these values).

Fig. 2
figure 2

a An analytical temperature field calculated by Eq. (1) and b its fusion zone

Full size image

Finite Element Thermal Analysis

A finite element thermal analysis was performed by ABAQUS (ABAQUS/CAE 2021; Dassault Systems Simulia Corp.) with self-developed subroutines on a desktop computer with an Intel Core i7 CPU and an NVIDIA GeForce RTX 3060. As shown in Fig. 3a, the model consisted of two parts, POWDER and SOLID, with different material properties. The section that the heat source passes through was divided into the 8-node linear heat transfer brick elements with a mesh size of 5.0 µm, while the rest of the model was composed of coarser elements.

Fig. 3
figure 3

a Model for the finite element thermal analysis, b a simulated temperature field, and c a two-dimensional array \({\text{FZ}}_{{{\text{sim}}}}\) with components 1 and 0 representing the fusion zone and the unmelt area, respectively

Full size image

The governing equation of heat transfer in the L-PBF process is as follows:

$$ \begin{array}{*{20}c} {\rho C_{p} \frac{\partial T}{{\partial t}} = \nabla \left( {k\nabla T} \right) + Q} \\ \end{array} $$
(2)

in which \(\rho \) is density, \({C}_{p}\) is specific heat capacity, \(k\) is thermal conductivity, \(T\) is temperature, and \(t\) is time. \(Q\) includes the heat losses by natural convection and radiation, \({q}_{{\text{c}}}\) and \({q}_{{\text{r}}},\) as well as the heat input by laser \({q}_{{\text{laser}}}\) with a heat source model, which is detailed in Sect. “Heat source models.” The center of the heat source model moved in the x direction through point O at \(\left( {x, \;y, \;z} \right) = \left( {0, \;0, \;30} \right)\) so that the plane at \(y=0.0\) was considered as the symmetric boundary. The heat losses by natural convection and radiation, \({q}_{{\text{c}}}\) and \({q}_{{\text{r}}},\) defined as the following equations are considered on the top surface:

$$ \begin{array}{*{20}c} {q_{c} = h_{{\text{c}}} \left( {T - T_{{{\text{amb}}}} } \right)} \\ \end{array} $$
(3)
$$ \begin{array}{*{20}c} {q_{r} = \sigma \varepsilon \left( {T^{4} - T_{{{\text{amb}}}}^{4} } \right)} \\ \end{array} $$
(4)

where \({h}_{{\text{c}}}\), \(\sigma \), \(\varepsilon \), and \({T}_{{\text{amb}}}\) are the convection coefficient, Stefan–Boltzmann constant, emissivity, and ambient temperature, respectively. The values of these physical constants are listed in Table 2. The other walls were set to be adiabatic, and the model width and depth were sufficiently large so that the effects of such boundary conditions were negligible during the solidification process.

Table 2 Material properties of IN738LC and physical constants
Full size table

The material properties of IN738LC were extracted from the simulation software package JMatPro (Sente Software Ltd.), which calculates a wide range of material properties. The constant values of \(\rho \), \(L\), \({T}_{{\text{s}}}\) and \({T}_{{\text{l}}}\) are also listed in Table 2, and the temperature-dependent specific heat capacity, \({C}_{p}\), and thermal conductivity, \(k\), are plotted in Fig. 4. As in the previous study [35], the latent heat \(L\) associated with phase transformation was modeled as a normal distributed equivalent heat capacity from \({T}_{{\text{s}}}\) to \({T}_{{\text{l}}}\) (see Fig. 4). As for \(k\), the value for the liquid phase (473 W/mK) was set to be more than 10 times larger than that for the solid phase. This can be understood by referring to the effective thermal conductivity for another nickel superalloy, Inconel 718, reported by Ladani et al. [26]. According to their study, the effective thermal conductivity was determined by considering the Marangoni effect (the surface tension gradient-driven flow) in the melting pool and validated by the comparison between the finite element modeling and experimental measurements. In addition, the thermal conductivity of the powder was initially set to be 5% of the bulk counterpart (= \(0.05 k\)), and gradually approached that of SOLID as the temperature rose from \({T}_{{\text{s}}}\) to \({T}_{{\text{l}}}\) (see the dotted light gray line in Fig. 4). After the element temperature in the POWDER part exceeded \({T}_{{\text{l}}}\), the material properties of the element were changed to those of the SOLID part. The density of the POWDER part was \(0.60 \rho \) even after the solidification to maintain the mass balance. The sensitivities of \(\kappa \) in \({k}_{{\text{POWDER}}}=\kappa {k}_{{\text{SOLID}}}\) and \(r\) in \({\rho }_{{\text{POWDER}}}=r{\rho }_{{\text{SOLID}}}\) were studied in the range of 0.01–0.10 and 0.50–0.65, respectively. The simulation error due to the changes in \(\kappa \) and \(r\) was found to be less than 3%, confirming that its effect on the simulated melt pool was limited.

Fig. 4
figure 4

Temperature-dependent heat capacity and thermal conductivity of IN738LC. The light gray dotted line indicates the thermal conductivity of powder, and the red and blue broken lines are the solidus and liquidus temperatures, respectively

Full size image

As shown in Fig. 3b, a transient temperature field was obtained through the thermal analysis. In addition, a state variable with an initial value of 0 was assigned to each integral point, and the value was changed to 1 when the temperature exceeded \({T}_{{\text{l}}}\) through the thermal analysis. The state variables at \(x=0\) were extracted as the cross-sectional array and referred to as \({\text{FZ}}_{{{\text{sim}}}}\) (see Fig. 3c).

Heat Source Models

The volumetric heat flux distributions by the four models are illustrated in Fig. 5. The conical model is represented as the following equation [24], and its heat flux distribution is illustrated in Fig. 5a:

$$ \begin{array}{*{20}c} {Q\left( {x, y, z} \right) = \frac{6\eta P}{{\pi h(r_{e}^{2} + r_{e} r_{i} + r_{i}^{2} )}}{\text{exp}}\left( { - 2\frac{{x^{2} + y^{2} }}{{r_{0}^{2} }}} \right)} \\ \end{array} $$
(5)
$$ \begin{array}{*{20}c} {r_{0} = r_{e} + \frac{z}{h}\left( {r_{e} - r_{i} } \right)} \\ \end{array} $$
(6)

where \({r}_{e}\) and \({r}_{i}\) are the radius at the top and bottom, respectively, and \(h\) corresponds to the height of the conical shape. In this study, the constraint \({r}_{e}\ge {r}_{i}\) was set. When \({r}_{e}={r}_{i}\), Eqs. (5) and (6) represent the cylinder model. Applying depth-decaying functions to the cylinder model yields the linear decaying model in Eq. (7) [26] and the exponentially decaying model in Eq. (8) [7, 27, 28]:

$$ \begin{array}{*{20}c} {Q\left( {x,\; y,\; z} \right) = \frac{2\eta P}{{\pi r_{e}^{2} }}{\text{exp}}\left( { - 2\frac{{x^{2} + y^{2} }}{{r_{e}^{2} }}} \right) \cdot \frac{2}{{h^{\prime}}}\left( {1 - \frac{z}{{h^{\prime}}}} \right)} \\ \end{array} $$
(7)
$$ \begin{array}{*{20}c} {Q\left( {x,\; y,\; z} \right) = \frac{2\eta P}{{\pi r_{e}^{2} }}{\text{exp}}\left( { - 2\frac{{x^{2} + y^{2} }}{{r_{e}^{2} }}} \right) \cdot \frac{1}{h\prime \prime }\exp \left( { - \frac{z}{h\prime \prime }} \right)} \\ \end{array} $$
(8)

where both \(h{\prime}\) and \(h{\prime}{\prime}\) are the shape parameters to define the decaying depth. Figure 5b and c shows the heat flux distributions of the linearly and exponentially decaying models, respectively.

Fig. 5
figure 5

Volumetric heat flux distributions by different source models: a conical, b linearly decaying, c exponentially decaying, and d double ellipsoidal (\({r}_{e}=30 \mu {\text{m}}\), \({r}_{i}=20 \mu {\text{m}}\) for the conical model, \({r}_{e}=50 \mu {\text{m}}\) for the linearly and exponentially decaying models, \(h = h^{\prime } = h^{\prime \prime } = 75 \;{\mu m},\;a_{f} = 30\; {\mu m},\;a_{r} = 100 \;{\mu m},\;b = 50\; {\mu m},\;\) and \(c = 75\; {\mu m}\))

Full size image

As described in Sect. “Introduction”, the most well-known heat source model, often used in the L-PBF process analysis, is the double ellipsoidal model proposed by Goldak et al. [19]. As shown in Fig. 5d, this model is composed of two ellipsoids with different parameters forward and backward with respect to the moving direction (\(x\)-direction in the figure):

$$ \begin{array}{*{20}c} {Q_{f} \left( {x,\; y,\; z} \right) = \frac{6\sqrt 3 \eta P}{{\pi \sqrt \pi a_{f} bc}}f_{f} \exp \left( { - 3\left( {\frac{{x^{2} }}{{a_{f}^{2} }} + \frac{{y^{2} }}{{b^{2} }} + \frac{{z^{2} }}{{c^{2} }}} \right)} \right)} \\ \end{array} $$
(9)
$$ \begin{array}{*{20}c} {Q_{r} \left( {x, \;y,\; z} \right) = \frac{6\sqrt 3 \eta P}{{\pi \sqrt \pi a_{r} bc}}f_{r} \exp \left( { - 3\left( {\frac{{x^{2} }}{{a_{r}^{2} }} + \frac{{y^{2} }}{{b^{2} }} + \frac{{z^{2} }}{{c^{2} }}} \right)} \right)} \\ \end{array} $$
(10)
$$ \begin{array}{*{20}c} {f_{f} + f_{r} = 2} \\ \end{array} $$
(11)

where \({a}_{f}\) or \({a}_{r}\), \(b\), and \(c\) are equivalent to the radius in the \(x\)-, \(y\)-, and \(z\)-axis directions, respectively. Based on previous experimental and simulation experiences [19, 22, 23, 35], the constraint \({a}_{f}\le {a}_{r}\) was placed.

Note that the movement of the heat source with scanning velocity \(v\) in the \(x\) direction is just a coordinate transformation, which can be expressed by replacing \(x\) in the above heat source equations with \(\left( {x - vt} \right)\).

Shape Parameter Calibration

The shape parameters were calibrated using the procedure developed in our previous study [35]. As described in Secs. “Bead-on-plate test” and “Finite element thermal analysis,” the arrays \({\text{FZ}}_{{{\text{obs}}}}\) and \({\text{FZ}}_{{{\text{sim}}}}\) were prepared by the experimental bead-on-plate test and the thermal analysis, respectively. As shown in Figs. 1a and 3c, its components 1 and 0 represent the fusion zone and the unmelted area, respectively. Thus, by taking the component difference \({\text{FZ}}_{{{\text{sim}}}}\) and \({\text{FZ}}_{{{\text{obs}}}}\), the simulation error \(\varepsilon\) was defined as follows:

$$ \begin{array}{*{20}c} {\varepsilon = 100\frac{{\mathop \sum \nolimits_{i = 1} \mathop \sum \nolimits_{j = 1} \left| {{\text{FZ}}_{{{\text{sim}}, \;ij}} - {\text{FZ}}_{{{\text{obs}},\;ij}} } \right|}}{{\mathop \sum \nolimits_{i = 1} \mathop \sum \nolimits_{j = 1} \left| {{\text{FZ}}_{{{\text{obs}},\;ij}} } \right|}}} \\ \end{array} $$
(12)

where \(i\) and \(j\) are the indices in the \(y\) and \(z\) directions in the arrays. Since \({\text{FZ}}_{{{\text{sim}}}}\) depends on the heat source models with the shape parameters, the optimal shape parameter values minimize \(\varepsilon \). Thus, as shown in Fig. 6, the optimal values were searched through 30 iterations of the thermal analysis with different values of the shape parameters proposed by Bayesian optimization. As in our previous study [35], the automatic relevance determination (ARD) Matérn 5/2 kernel and expected improvement (EI) were used as kernel and acquisition functions in Bayesian optimization, respectively. The procedure was executed by using Python 3.8 with the Bayesian optimization module [38].

Fig. 6
figure 6

Calibration procedure of the shape parameters

Full size image

Results

Bead-on-Plate Test

Figure 7 shows the cross-sectional SEM images with different laser powers and scanning velocities. The boundary of the fusion zone can be seen by the difference of microstructure. Obviously, the depth of the fusion zone was increased with the increase in the laser power (see Fig. 7a, b, and c) and the decrease in the scanning velocity (see Fig. 7f, e, and d). As indicated green arrows in Fig. 7b, c, e, and f, the dents were found at the base of the bead under the conditions of \(P>400\) W and \(v>1500\) mm/s. This is often observed in welding processes at fast scanning velocities, and is referred to as humping [39]. Also, as indicated by the white arrow in Fig. 7f, a large porosity was formed at the bottom of the fusion zone. As shown in Fig. 8, the depth \({d}_{{\text{obs}}}\) and area \({A}_{{\text{obs}}}\) of the fusion zone increased almost linearly with \(P/v\), and the width \({w}_{{\text{obs}}}\) also tended to increase with the increase of \(P/v\).

Fig. 7
figure 7

Cross-sectional SEM images for the IN738LC beads at different \(P\) and \(v\): a 300 W, 1500 mm/s, b 500 W, 1500 mm/s, c 700 W, 1500 mm/s, d 400 W, 1000 mm/s, e 400 W, 1500 mm/s, and f 400 W, 2000 mm/s. The yellow broken line in (a) indicates the outline of the fusion zone. The green and white arrow indicate the dents at the base of a bead and a large porosity at the bottom of the fusion zone, respectively

Full size image
Fig. 8
figure 8

a Depth \({d}_{{\text{obs}}}\), b width \({w}_{{\text{obs}}}\), and c area \({A}_{{\text{obs}}}\) of the fusion zone in the bead-on-plate test plotted with \(P/v\). The color and size of the plotted points indicate the laser power and scanning velocity, respectively

Full size image

Estimation of Effective Absorptivity

The fusion zone for each laser scanning condition was calculated using the analytical solution of Eq. (1). Figure 9a shows the fusion zone with \(\eta \) of 0.1, 0.5, and 0.9 at \(P=\) 300 W and \(v=\) 1500 mm/s. It is obvious that as \(\eta \) increased, the fusion zone was larger due to the higher input energy. Thus, the bisection method was used to determine the effective absorptivity \({\eta }^{{\text{BS}}}\) where the area \({A}_{{\text{as}}}\) is consistent with \({A}_{{\text{obs}}}\), or \(\left|F\left(\eta \right)\right|<0.005\). As shown in Fig. 9b, \({\eta }^{{\text{BS}}}\) for \(P=\) 300 W and \(v=\) 1500 mm/s was found to be 0.65.

Fig. 9
figure 9

a Cross-sectional fusion zones in the case of \(\eta =\) 0.1, 0.5, and 0.9 at \(P=\) 300 W and \(v=\) 1500 mm/s. b \(F\left(\eta \right)\) in the process of determining \({\eta }^{{\text{BS}}}\) by the bisection method

Full size image

In the same way, \({\eta }^{{\text{BS}}}\) was determined for all the conditions. As shown in Fig. 10a, \({\eta }^{{\text{BS}}}\) tended to increase as the laser energy density \(P/v\) increased. As in the previous studies [14, 31], this result is reasonable because the more the laser energy density increases, the more concave the melting pool surface becomes, and consequently the more multiple reflections of the laser beam occur. Unfortunately, to the best of the authors’ knowledge, no previous studies have reported the effective absorptivity of IN738LC powder measured by calorimetric measurement. However, according to the study using calorimetric measurement by Boley et al. et al. [33], the absorptivity of nickel-based alloy IN625 was almost the same as that of stainless steel 316L (SS316L). In addition, Trapp et al. [31] measured the absorptivity of SS316L with different laser scanning conditions using calorimetric measurement. Although the composition differs between IN738LC and SS316L, their results are also shown in Fig. 10a. This result shows that the effective absorptivity estimated from Eq. (1) is close to that measured by calorimetric measurement.

Fig. 10
figure 10

a \({\eta }^{{\text{BS}}}\) against \(P/v\) and b \({\eta }^{{\text{MLR}}}\) against \({\eta }^{{\text{BS}}}\). The color and size of the plotted points indicate the laser power and scanning velocity, respectively. The gray square plots in (a) were the absorptivity measured by Trapp et al. [31]

Full size image

To represent \({\eta }^{{\text{BS}}}\) as a function of the laser scanning conditions, multiple linear regression (MLR) analysis was conducted with \(P\), \(v\), \(P/v\), and \(v/P\) as the explanatory variables and \({\eta }^{{\text{BS}}}\) as the objective variable. The appropriate explanatory variables were selected with the leave-one-out cross-validation (LOOCV) method [40] to improve generalization performance of the regression model by avoiding over-fitting. The derived regression model is as follows:

$$ \begin{array}{*{20}c} {\eta^{{{\text{MLR}}}} = 9.39 \times 10^{ - 5} v + 1.01\frac{P}{v} + 2.71 \times 10^{ - 2} \frac{v}{P} + 9.94 \times 10^{ - 2} } \\ \end{array} $$
(13)

As shown in Fig. 10b, where \(\eta^{{{\text{MLR}}}}\) is plotted against \(\eta^{{{\text{BS}}}}\), the plots line up near the solid black line of \({\eta }^{{\text{MLR}}}={\eta }^{{\text{BS}}}\). The mean absolute error between \({\eta }^{{\text{BS}}}\) and \({\eta }^{{\text{MLR}}}\) was 0.030.

Screening Heat Source Models

As an example, the calibration results of the shape parameters for the conical model at \(P=400\) W and \(v=1000\) mm/s are summarized in Fig. 11 (see also the online supplemental movies, S1, S2, S3, S4 and S5). Figure 11a and b shows the shape parameters \(h\) and \({r}_{e}\) against \({r}_{i}\) for each trial, and the plots are connected by black straight lines in the order of the trials. Each time new values were determined for the shape parameters by the Bayesian optimization, the thermal analysis was executed using the heat source model (Fig. 11c) to obtain \({{\text{FZ}}}_{{\text{sim}}}\) (Fig. 11d). Then, the simulation error \(\varepsilon \) was determined and plotted with the trial number as a red line in Fig. 11e. In addition, the blue line in the figure shows the minimum \(\varepsilon \) up to the latest trial. The plot color in Fig. 11a and b also corresponds to the value of \(\varepsilon \). As shown in Fig. 11a and b, the parameter search was mainly focused on the area with small errors while also examining the entire region due to the balance between exploitation and exploration of EI in the Bayesian optimization. Consequently, as indicated by the intersection of the red lines in Fig. 11a and b and by the dotted broken line in Fig. 11e, the best shape parameters to minimize \(\varepsilon \) were determined within 30 trials. The fusion zone \({{\text{FZ}}}_{{\text{sim}}}\) by the thermal analysis with the best shape parameters is shown in Fig. 11d and was in good agreement with \({{\text{FZ}}}_{{\text{obs}}}\) by the SEM observation.

Fig. 11
figure 11

Shape parameter calibration for the conical model at \(P=400\) W and \(v=1000\) mm/. a Shape parameter \(h\) against \({r}_{i}\), and b \({r}_{e}\) against \({r}_{i}\) for each trial. The point of intersection among the red lines indicates the best shape parameters. c Heat flux distribution of the conical model with the best parameters. d \({FZ}_{obs}\), \({FZ}_{sim}\), and the difference between them. e Error with the number of trials. Red lines represent the error for each trial, and the blue line indicates the lowest error up to the latest trial (see also the online supplemental movies, S1, S2, S3, S4, and S5)

Full size image

In the same way, the shape parameters for the four heat source models were also calibrated for the representative five laser scanning conditions. The selected conditions were the 1st, 5th, 9th, 13th, and 17th in ascending order of \({A}_{{\text{obs}}}\) in the bead-on-plate test in Sect. “Bead-on-plate test.” Figure 12 shows a heat map representing the minimum \(\varepsilon \) for each model and condition after the calibration. The conical model had the smallest \(\varepsilon \) among the four heat source models in the four conditions except for the case of \(P=\) 400 W and \(v=\) 2000 mm/s. Consequently, as shown in the bottom row, labeled “Average for all conditions,” the exponentially decaying model tended to have the largest \(\varepsilon \), followed by the double ellipsoidal and linearly decaying models, and the conical model had the smallest.

Fig. 12
figure 12

Heat map of the minimum \(\varepsilon \) with the calibrated shape parameters each heat source model and laser scanning condition. The bottom row shows the mean \(\varepsilon \) for the five laser scanning conditions

Full size image

Figure 13 shows the temperature distributions at \(P=\) 400 W and \(v=\) 1000 mm/s simulated by the thermal analysis with the conical and double ellipsoidal models. The center of each heat source model passed through the \(yz\) plane including point O at \(t=0.20\) ms. The yellow dotted line in the figure represents the outline of the fusion zone from the bead-on-plate test in Sect. “Bead-on-plate test.” In the case of the double ellipsoidal model in Fig. 13b, the area below the fusion outline was directly heated by the bottom of the heat source (see Fig. 5d). The same tendency was observed in the results by the linearly and exponentially models in the supplementary image S6. Unlike these depth-decaying models, the areas outside the fusion outline shown in Fig. 13a were barely heated by the conical model, which is a depth-constant model and does not attenuate in the depth direction (see Fig. 5a). Some CFD-based thermal simulations [12,13,14] showed that the melting surface was hollowed out by the laser so that the bottom of the melting pool was directly heated, but the region below the heated surface was not directly heated because the laser beam is impermeable to metal. Thus, such an area below the fusion outline should be heated indirectly and gradually by thermal diffusion after the heat source has completely passed through. These results conclude that the depth-constant conical heat source model is the best approximation of laser melting in the L-PBF process than the depth-decaying models.

Fig. 13
figure 13

Temperature distributions at \(P=400\) W and \(v=1000\) mm/s simulated by the thermal analysis with a the conical and b double ellipsoidal heat source models with the calibrated shape parameters. The center of the heat source model is on the plane at \(t=0.20\) ms. The yellow dotted lines represent the outline of the fusion zone from the bead-on-plate test. See also the online supplemental Fig. S6 for the results with the other heat source models

Full size image

Functionalization of Shape Parameters as Laser Scanning Conditions

Figure 14 shows the shape parameters for the conical model with all the laser scanning conditions. It is obvious that the height parameter \({h}^{{\text{BO}}}\) increases linearly with the increase of \(P/v\) (see Fig. 14a). On the other hand, the dependency of the radial parameters \({r}_{e}^{{\text{BO}}}\) and \({r}_{i}^{{\text{BO}}}\) on the laser scanning conditions was unclear in Fig. 14b and c. Thus, in the same way as \({\eta }^{{\text{MLR}}}\) derived as a function of \(P\) and \(v\) in Sect. “Estimation of effective absorptivity,” MLR analysis with LOOCV was conducted with \(P\), \(v\), \(P/v\), and \(v/P\) as the explanatory variables and each shape parameter for the conical model as the objective variable. As a result, the following regression equations were obtained:

$$ \begin{array}{*{20}c} {h^{{{\text{MLR}}}} = 0.0523P + 705\frac{P}{v} + 24.7\frac{v}{P} - 178} \\ \end{array} $$
(14)
$$ \begin{array}{*{20}c} {r_{e}^{{{\text{MLR}}}} = 0.0650P + 109\frac{P}{v} + 11.4\frac{v}{P} - 64.1} \\ \end{array} $$
(15)
$$ \begin{array}{*{20}c} {r_{i}^{{{\text{MLR}}}} = 0.143P - 0.0285v - 156\frac{P}{v} - 4.61\frac{v}{P} + 73.4} \\ \end{array} $$
(16)
Fig. 14
figure 14

Calibrated shape parameters plotted against \(P/v\): a \({h}^{{\text{BO}}}\), b \({r}_{e}^{{\text{BO}}}\), and c \({r}_{i}^{{\text{BO}}}\) for the conical model. The color and size of the plotted points indicate the laser power and scanning velocity, respectively

Full size image

These equations suggest that the shape parameters are not just proportional to the energy density \(P/v\), but represented as nonlinear equations of \(P\) and \(v\). As shown in Fig. 15a, b, and c, the shape parameters \({h}^{{\text{MLR}}}\), \({r}_{e}^{{\text{MLR}}}\), and \({r}_{i}^{{\text{MLR}}}\) were plotted against \({h}^{{\text{BO}}}\), \({r}_{e}^{{\text{BO}}}\), and \({r}_{i}^{{\text{BO}}}\), respectively. The plots lined up near the solid black line imply that the regression models. The mean absolute errors between \({h}^{{\text{BO}}}\) and \({h}^{{\text{MLR}}}\), \({r}_{e}^{{\text{BO}}}\) and \({r}_{e}^{{\text{MLR}}}\), and \({r}_{i}^{{\text{BO}}}\) and \({r}_{i}^{{\text{MLR}}}\) were 8.16, 7.45, and 5.50 µm, respectively.

Fig. 15
figure 15

Shape parameters as a function of \(P\) and \(v\) against the calibrated ones: a \({h}^{{\text{MLR}}}\) against \({h}^{{\text{BO}}}\), b \({r}_{e}^{{\text{MLR}}}\) against \({r}_{e}^{{\text{BO}}}\), and c \({r}_{i}^{{\text{MLR}}}\) against \({r}_{i}^{{\text{BO}}}\). The color and size of the plotted points indicate the laser power and scanning velocity, respectively

Full size image

Now that the conical heat source model was finally represented as laser scanning conditions, the thermal analysis with the MLR model was performed under each of all the 17 conditions in Fig. 2. The simulation error was evaluated based on Eq. (12), and the mean error averaged over all the 17 conditions was 16.6 \(\pm \) 9.0%. Also, comparing the fusion zones of the thermal analysis and the bead-on-plate test, the mean absolute errors of the depth and width of were 5.4 \(\pm \) 4.0 µm, 4.6 \(\pm \) 2.6 µm, respectively.

In order to confirm the validation of the MLR heat source model at arbitrary laser scanning conditions, the thermal analysis with the model was performed with the \(P=300\) W, \(v=700\) mm/s and \(P=300\) W, \(v=1100\) mm/s, which are not included in the 17 conditions used for the MLR analysis. In addition, assume that if there are no MLR heat source model and no re-calibration for the new conditions, it is practical to compromise by carrying out the thermal analysis by substituting the alternative shape parameter values that have been already calibrated for a close laser scanning condition. Thus, for comparison, the thermal analysis for each new condition was also tried by substituting the alternative shape parameter values calibrated for \(P=300\) W and \(v=1000\) mm/s. Figure 16a and b shows, from left to right, the fusion zone \({{\text{FZ}}}_{{\text{obs}}}\) from the bead-on-plate test, the fusion zone \({{\text{FZ}}}_{{\text{sim}}}\) and difference by the analysis substituting the alternative shape parameters instead of re-calibration, and the fusion zone \({{\text{FZ}}}_{{\text{sim}}}\) and difference from the analysis using the MLR model. The figure shows that \({{\text{FZ}}}_{{\text{sim}}}\) was more identical to \({{\text{FZ}}}_{{\text{obs}}}\) using the MLR model than substituting the alternative shape parameters. According to Eq. (12), the error using the MLR was 18.2% for \(P=300\) W, \(v=700\) mm/s, and 7.6% for \(P=300\) W, \(v=1100\) mm/s. On the other hand, the errors substituting the alternative parameters were 45.2% for \(P=300\) W, \(v=700\) mm/s and 18.1% for \(P=300\) W, \(v=1100\) mm/s, which were more than twice larger than those with the MLR model. As a result, without additional calibrations, the MLR heat source model enables the more accurate thermal analysis under arbitrary laser scanning conditions in the process window.

Fig. 16
figure 16

\({{\text{FZ}}}_{{\text{obs}}}\), \({{\text{FZ}}}_{{\text{sim}}}\) and their difference using the heat source model with the shape parameters calibrated for \(P=300\) W, \(v=1000\) mm/s, and those using the MLR model: a \(P=300\) W, \(v=700\) mm/s and b \(P=300\) W, \(v=1100\) mm/s

Full size image

Discussion

These results provide a heat source model that can be used at arbitrary \(P\) and \(v\) within the process window. One advantage of MLR over other regression methods (support vector machine, random forest algorithm, and so on) is that it can be expressed as a polynomial of explanatory variables. In other words, by just transcribing the polynomials of shape parameters (Eqs. (14)–(16)) into the moving heat source subroutine, the thermal analysis becomes valid within the process window. In addition, LOOCV prevented over-training on the training data so that Eqs. (14)–(16) can be valid in any laser scanning condition within the process window. It should be noted that the process window here is within the laser scanning conditions used as training data for the MLR analysis, i.e., interpolation between the conditions is valid, but extrapolation is not guaranteed. That is, the process window in which the MLR model is valid is specifically as follows:

$$ \begin{gathered} 100 \le P \le 700 \hfill \\ \begin{array}{*{20}c} {500 \le v \le 2000} \\ \end{array} \hfill \\ 0.133 < P/v < 0.5 \hfill \\ \end{gathered} $$
(17)

As described in Introduction, the effective absorptivity and shape parameters properly expressed as a function of the laser scanning conditions eliminate the need for re-calibration which limits the thermal analysis for an arbitrary laser scanning condition. In the following, a case study utilizing the MLR heat source model is presented to search for the optimal laser scanning condition that achieve the target values. One of the target parameters was the depth of the fusion zone, and the target depth \({d}^{{\text{Target}}}\) was set to 120.0 µm. The other target parameter was the cooling rate \({\text{d}}T/{\text{d}}t\) at \(T_{{\text{l}}}\). As shown in Fig. 17b, the cooling rate was largest at the bottom of the fusion zone and smaller at shallower depth. The target value \({{{CR}}}^{{\text{Target}}}\) at point P, 30 µm shallow from the bottom, was set at 1.1 \(\times \) 106 K/s. An indicator of how close the simulated values \(d\) and \({{CR}}\) is to these target values \(d^{{{\text{Target}}}}\) and \({{CR}}^{{{\text{Target}}}}\) was defined as follows:

$$ \begin{array}{*{20}c} {G = \frac{{\left| {d - d^{{{\text{Target}}}} } \right|}}{{d^{{{\text{Target}}}} }} + \frac{{\left| {{{CR}} - {{CR}}^{{{\text{Target}}}} } \right|}}{{{{CR}}^{{{\text{Target}}}} }}} \\ \end{array} $$
(18)
Fig. 17
figure 17

a Procedure to search the optimal laser scanning condition to meet the target values. b Cooling rate at \({T}_{{\text{l}}}\) at the \(yz\) plane. c Function \(G\), e depth \(d\), and g cooling rate \({{CR}}\) through the search. The contour maps for d \(d\) and f \({{CR}}\)

Full size image

As shown in Fig. 17a, the optimal laser scanning condition \(P\) and \(v\) to minimize \(G\) was searched within the process window in Eq. (17) by the iterative thermal analysis with the Bayesian optimization using the same kernel and acquisition functions as the shape parameter calibration. Figure 17c, e, and g shows \(G\), d, and \({{CR}}\) against the number of searches, respectively. As indicated by the black dotted lines in the figure, the 42nd thermal analysis out of 44 searches resulted in a minimum value of 0.0233 for \(G\), and the d and \({{CR}}\) at 42nd were respective 122.5 µm and 1.097 \(\times \) 106 K/s, and significantly close to the target values indicated by the horizontal broken lines in Fig. 17e and g). Figure 17d and f shows the contour map of d and \({{CR}}\), respectively. As indicated by the white circles in the figure, the laser scanning conditions \(P\) and \(v\) suggested by the Bayesian optimization distribute around the target values as well as around the edge of the process window. The red plot in the figure represents the 42nd condition, which suggests the optimal conditions \({P}^{{\text{opt}}}=\) 423 W and \({v}^{{\text{opt}}}=\) 1348 mm/s to achieve the target values. In terms of directions for future research, it would be promising to couple this thermal analysis with other analyses (e.g., thermo-mechanical analysis) and execute process optimization with another indicator (such as stress distribution).

Conclusions

The current study set out to develop a heat source model for finite element thermal analysis that can be applied to a wide process window of L-PBF. A secondary objective was to fairly and quantitatively evaluate four heat source models proposed so far and find the best one for the L-PBF process. Thanks to a combination of the quantification of the error and the efficient calibration procedure using the Bayesian optimization, it was found that the conical model is most suitable for the thermal analysis of L-PBF rather than depth-decaying models including the double ellipsoidal one. This approach will prove useful in assessing other heat source models not addressed in this study. The MLR analysis with LOOCV successfully represented the effective absorptivity and shape parameters for the conical model as a function of the laser power and scanning velocity. The MLR heat source model allowed the thermal analysis to be valid and precise under any laser scanning condition within the process window. While the mean simulation error using the MLR model was 16.6 \(\pm \) 9.0%, the mean absolute errors of the fusion depth and width were 5.4 \(\pm \) 4.0 µm and 4.6 \(\pm \) 2.6 µm, respectively. The need for re-calibration for each arbitrary condition has been removed by the MLR heat source model, which enables process optimization as demonstrated in the case study.