Introduction

Additive manufacturing (AM) is a revolutionary manufacturing process with clear advantages such as lower energy usage, minimum scrap waste, lower buy-to-fly ratio and shorter lead time to market.[1] The laser powder-bed fusion additive manufacturing (LPBF) process is one of the most successful process to realize lightweight and cost-effective production of complex, high-performance end-use parts. The widespread and economical use of metal AM relies on the ability to predict and control microstructures and resulting mechanical properties.[2] To date, the cost and time associated with process development for LPBF of new alloys is high and deployed processes were not optimized. A large gap exists in the AM process development to link process parameters to microstructure and final part performance. This prevents the full exploration of the component design space to the realization of full benefits of AM components.

Process maps for laser-engineered net shaping (LENS) process were developed.[3,4,5] However, the extension to LPBF of those process maps is constrained by computational resources and consolidated solidification theory with insufficient validation data, due to much higher spatial and temporal resolution specific to LPBF. By varying particle size distribution (PSD) of powder, laser parameters (power, scan speed, beam diameter, and hatch spacing), metal alloys can experience large range of solidification conditions,[6] namely temperature gradient (G) and solidification velocity (V), from ~ 103 to 106 K/m and ~ 0.001 to 10 m/s, respectively. The resultant solidification microstructure will vary dramatically based on the process parameters.[7] Current approaches to predict thermal history in LPBF are solving mainly for the energy equation using analytical approaches, e.g., Rosenthal solution[6] or commercially available software solving the energy equation without the fluid flow.[8] Recently, open-source software[9,10,11] and few commercially available software[12,13] that include fluid dynamics effects were used for LPBF process simulations.

However, most of the studies on predicting grain morphology based on solidification map correlations were conducted for heat-transfer-solidification-only (HTS) models, i.e., without including the fluid dynamics effects. This was partly due to either the lack of comprehensive solidification models that include the fluid dynamic effects or excessive computational time required to conduct coupled fluid dynamics and solidification simulations (HTS-FD).[12,13] To drastically reduce the simulation time and increase the use of fluid dynamics-based models in industry, phase-change models must be implemented in massively parallel CFD codes to simulate the LPBF process as a typical LPBF process involves ~30 μm powder layer thickness with melt-pool width ~ 150 μm width moving at speeds as high as a few meters per second; while melting-solidification phenomena occur at sub-micron scale. Most massively parallel CFD commercial codes, which were developed for general applications, were not used for AM simulations. Recently, several codes were developed for AM simulation.[9,10,11,14,15] However, their availability is limited due to software costs or are used internally at research centers. Truchas, is an open-source massively parallel CFD software that was developed under the Advanced Simulation and Computing (ASC) program at Los Alamos National Laboratory (LANL) for the simulation of metal-casting processes.[16] However, Truchas has been used to simulate the heat transfer and solidification of electron beam additive manufacturing process[17,18,19,20] without including fluid dynamics effects.

Thermo-capillary forces due to the surface tension variation and large temperature gradients across the molten metal surface, are the main driving force for fluid flow during LPBF. To simplify the modeling of surface tension effects on free-deformable surfaces, the flat-surface assumption for thermo-capillary forces was extensively used for welding, including those laser-based.[21,22,23,24,25,26] In flat-surface models, the molten alloy surface was held fixed by imposing no-flow boundary condition in the direction normal to the free surface. As reported,[26] most of the flat-surface laminar models used for welding employed multiplier factors for the thermal conductivity and the viscosity to obtain a good agreement with the measured weld pool shapes.[23,24,27,28]

The properties of a part produced by LPBF depend strongly on the properties of each single track and each single layer.[29] Considering the track-by-track and layer-by-layer deposition, the microstructure of the LPBF was found to exhibit a layer-wise layout within a vertical cross-section through the build, in which apparent molten pool outline curves can be easily differentiated.[30] Thus, the analysis of the microstructure distribution within a single track is a focus of this study.

The microstructure types and its length scales were shown to be governed by the following variables[7]: local temperature gradient (G), solidification velocity (V), and combinations thereof: the cooling rate (dT/dt, or GV) and G/V. While G/V can be used to assess the type of the solidification front (planar, cellular, or dendritic), GV (cooling rate) controls the size of the solidification structure.[17,18,31,32,33,34] Although the laser-based AM processes feature nominally a smaller beam size, e.g., 100 μm, the solidification-related variables G and V were found to vary within the solidified melt-pool, which is of sub-millimeter length-scale. From heat-transfer and solidification simulations that excluded fluid dynamics effects, it was found that both G and V increase as the melt-pool becomes smaller with increasing the laser scan speed, U, for nickel-base superalloys.[35] On the other hand, G decreases as the melt-pool becomes larger with increasing laser power, P, for a constant laser scan speed, U. It was found that low V’s and high G’s exist at the bottom of the melt-pool, whereas a higher V’s and lower G’s exist close to the melt-pool surface.[34] To date, the effect of the fluid flow on the spatial distribution of the V and G variables has yet to be fully investigated.

In this work, a scalable open source code, Truchas, was used. As a step toward a comprehensive multi-physics LPBF model, single-track laser fusion (STLF) simulations were conducted for the nickel-base superalloy (IN625) assuming a flat-free surface of the molten alloy. Although this assumption is expected to underpredict the meltpool depth for IN625 (which has negative surface tension coefficient),[26] important insights can be obtained on the effect of fluid flow on the solidification during LPBF. The variation of the microstructure-related variables, G, V, G/V, and GV over the solidified melt-pool was assessed, to understand the microstructure variation in LPBF process.

Constitutive Equations for Laser Fusion Model

The energy transport model implemented in Truchas accounts for alloy solidification. The incident heat flux, \( q^{\prime\prime}_{L} \), at any location (x, y) on the top surface from the laser beam is given by the Gaussian profile, as:

$$ q^{\prime\prime}_{\text{L}} \left( {x, y} \right) = \frac{{2 P_{\text{L}} }}{{\pi R^{2} }}e^{{ - 2\frac{{\left( {x - x_{o} \left( t \right)} \right)^{2} + \left( {y - y_{\text{o}} \left( t \right)} \right)^{2} }}{{R^{2} }}}} , $$
(1)

where PL is the laser power, R is the laser beam radius, and (xo(t), yo(t)) indicates the position of the center of the laser beam. Assuming an emission from a gray body, the total heat flux into the top surface—including losses due to natural/forced convection, thermal radiation, and evaporation, is given as:

$$ q^{\prime\prime} = A_{\lambda } q^{\prime\prime}_{\text{L}} - h_{\text{C}} \left( {T - T_{\text{A}} } \right) - \varepsilon_{\infty } \sigma \left( {T^{4} - T_{\text{A}}^{4} } \right) - q_{\text{evap}} , $$
(2)

where, \( A_{\lambda } \), \( \varepsilon_{\infty } \), hC, T, TA, σ, are the absorptivity of the alloy surface at the wavelength of the laser, total hemispherical emissivity of the top surface, heat transfer coefficient, surface temperature, ambient temperature, and Stefan–Boltzmann constant, respectively. The evaporation model was based on evaluating: (a) the saturated vapor pressure of liquid alloy using the Clausius–Clapeyron equation, and (b) the heat flux loss due to evaporation, qevap, of a mass flux given by the Hertz–Knudsen–Langmuir equation[36] as a function of the latent heat of evaporation, Levap, molecular weight, M, surface temperature, Ts, saturated vapor pressure, and universal gas constant, R, as:

$$ q_{\text{evap}} = L_{\text{evap}} \beta p_{\text{o}} { \exp }\left[ { - \left( {{\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {T_{\text{s}} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${T_{\text{s}} }$}} - {\raise0.7ex\hbox{$1$} \!\mathord{\left/ {\vphantom {1 {T_{\text{o}} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${T_{\text{o}} }$}}} \right)L_{\text{evap}} M/R} \right]\sqrt {\frac{M}{{2 \pi R T_{\text{s}} }}} , $$
(3)

where po is the reference ambient pressure, and saturation temperature, To = Tsat(po). Here, β is an empirical coefficient which was introduced to account for condensation effects.

A modified projection (or, fractional-step) algorithm was used to solve the Navier–Stokes.[37,38,39,40] To be concise, only the equations related to the implementation of the surface tension algorithm are given in this study. The surface tension force term in the tangential direction, FS(T), appears in the predictor step of the projection-based algorithm,[41] in which an intermediate velocity, u*, is computed from the momentum equation, as:

$$ \rho_{\text{o}} \left( {\frac{{{\mathbf{u}}^{*} - {\mathbf{u}}^{n} }}{\Delta t} + {\mathbf{u}} \cdot \nabla {\mathbf{u}}} \right) = \mu \nabla^{2} {\mathbf{u}}^{n} + \left( {\rho - \rho_{\text{o}} } \right){\mathbf{g}} - \rho_{\text{S}} {\mathbf{u}}^{n} \frac{{\partial g_{\text{L}} }}{\partial t} + {\mathbf{FS}}^{\left( T \right)} , $$
(4)

where \( \rho_{\text{o}} \) is the liquid alloy density, \( \Delta t \) is a time step, μ is the dynamic viscosity of the fluid, gL is the volumetric fraction of the liquid,[40] and superscripts “n” indicates the variables at the previous time level. All the properties are listed in the Appendix B. The fluid flow is considered to be in the laminar regime based on the Re number of approximately 780 (Appendix B). The Navier–Stokes equation is solved for all the computational cells in which the liquid fraction is above a given coherency threshold value, \( g_{\text{L}}^{\text{coh}} \):

$$ g_{\text{L}} \ge g_{\text{L}}^{\text{coh}} \;\left( {{\text{for}}\;{\text{fluid}}\;{\text{cells}}} \right) $$
(5)

The default value for \( g_{\text{L}}^{\text{coh}} \) in Truchas was considered to be 0.01. All the fluid flow simulations were conducted using this default value for the \( g_{\text{L}}^{\text{coh}} \), unless otherwise noted.

Surface tension varies linearly with temperature. The density variation with temperature is considered by using the Boussinesq approximation. The implementation in Truchas of the thermo-capillary model was validated[42] against analytical results of Reference 43 for a “rigid-lid” cavity problem. For a flat-surface normal to the Z-axis of the system of coordinate, the components of the FS(T) are non-zero only in the cells on the top surface and are given as:

$$ {\mathbf{FS}}^{\left( T \right)} = \frac{d\sigma }{dT} \left( {\frac{\partial T}{\partial x}, \frac{\partial T}{\partial y},0 } \right). $$
(6)

Because the interface is flat and constant (having zero curvature), the surface tension force has only a tangential component. Thus, the projection step is comprised of the following update of the total pressure, \( P = p + \rho \;g\;h \), and velocity:

$$ \nabla \cdot \left( {\sigma_{\text{p}} \nabla P^{n + 1} } \right) = \frac{1}{\Delta t}\nabla \cdot \left( {\rho g_{\text{L}} {\mathbf{u}}^{*} } \right)\;{\text{and}}\;\rho g_{\text{L}} \frac{{{\mathbf{u}}^{n + 1} - {\mathbf{u}}^{*} }}{\Delta t} = \sigma_{\text{p}} \nabla P^{n + 1} , $$
(7)

where the projection variable, \( \sigma_{\text{p}} \), is given by:

$$ \sigma_{\text{p}} = \frac{1}{\rho } . $$
(8)

Evaluation of Microstructure from LPBF Simulations

Solidification in LPBF process governs the size and morphology of microstructure features, such as primary arm spacing and the extent of micro-segregation, which in turn affect the properties of the additively manufactured parts. The microstructure types and its length scales were shown to be governed by the following variables[7]: local temperature gradient (G), solidification velocity (V), and combinations thereof: the cooling rate (dT/dt, or GV) and G/V. While G/V can be used to assess the type of the solidification front (planar, cellular, or dendritic), GV (cooling rate) controls the size of the solidification structure.[17,31,33,35] The temperature gradient (G) is calculated, as:

$$ G = \sqrt {G_{x}^{2} + G_{y}^{2} + G_{z}^{2} } , $$
(9)

where Gx, Gy and Gz are temperature gradients along X, Y and Z directions, respectively. The solidification velocity, or growth rate, is calculated as the ratio between cooling rate, dT/dt, and temperature gradient, as:

$$ V = \frac{dT}{dt}\frac{1}{G}. $$
(10)

In this study, the solidification microstructure-related variables (G, V, G/V, and G V) were evaluated at a fraction solid fS = 0.05.[44]

Correlations for primary dendrite arm spacing (PDAS)

PDAS has been an important measure to describe the microstructure characteristics. The dependence of PDAS, denoted by \( \lambda_{1} \), on G and V was reviewed by Ghosh et al.[34] for nickel-based superalloys, where experimental data was referenced.[45,46] The power law dependence based on Hunt[47] and Kurz and Fisher’s[48] models is the most widely applied solution for PDAS prediction for AM and it can be written in a general form, as

$$ \lambda_{1} = A_{\text{o}} G^{ - m} V^{ - n} , $$
(11)

where m and n are positive exponents. Analytical expressions for constant \( A_{\text{o}} \) as a function of material properties were derived by assuming spherical or ellipsoidal dendrite tips.[47,48] We have to keep in mind that in all of the PDAS models, such as the ones by Hunt and Kurz–Fisher models, the shape of the dendrite tips under severe solidification conditions and fluid convection may differ than the assumed shape used to derive these models. In Hunt’s model, \( A_{\text{o}} = A \left( {k\Gamma \Delta T_{\text{o}} D_{\text{L}} } \right)^{0.25} \), while in Kurz–Fisher’s model, \( A_{\text{o}} = 4.3 \left( {\frac{{\Gamma \Delta T_{\text{o}} D_{\text{L}} }}{k}} \right)^{0.25} \)where k is the equilibrium partition coefficient, \( \Gamma \) the Gibbs–Thomson coefficient, \( \Delta T_{\text{o}} \) the equilibrium freezing range, DL the diffusivity of the liquid. In these analytical models, the exponents in Eq. [11] were m = 0.5 and n = 0.25 for both Hunt and Kurz–Fisher models. Another power law for the PDAS dependence was shown to be simply \( \left( {G V} \right)^{ - m} \).[46]

Aside from these simple power law expressions of PDAS at high-cooling rates, Trivedi and Kurz[46] provided a more complex formulation. Trivedi and Kurz[49] summarized a dendritic growth model that accounts for attachment kinetics and composition variations at the interface at high growth rate, which accommodates the non-equilibrium effects during rapid solidification (Appendix C). Based on the conservation of dendrite tips, a dendrite tip radius selection criterion, a marginal stability criterion was adopted to resolve the unique tip radius and velocity under steady-state. The application of this model has been used to predict tip velocity, tip radius and solute trapping at high Peclet number conditions.[50] However, the Trivedi–Kurz (TK) model for PDAS has not been used for AM modeling; in part due to its much more complicated, non-linear analytical constitutive formulations.

Keller et al.[51] showed for IN625 that Hunt’s model yields PDAS values closer to the experimental values than Kurz–Fisher model, using G and V values from heat-transfer-and-solidification simulations, similar to the HTS model in this study. Since the differences between Hunt’s and Kurz–Fisher’s models are limited to the constants to the power laws, the general form of \( \lambda_{1} \left( {G, V} \right) = A_{\text{o}} G^{ - 0.5} V^{ - 0.25} \) is considered in this paper. The actual PDA estimates depend on the selection of material thermo-physical properties; some of them not readily available for multicomponent alloys as they are difficult to estimate from first principles and difficult to measure. Thus, where appropriate, some of the parameters in the PDAS models can be calibrated based on experimental data. For example, Ghosh et al.[35] used (GHTS, VHTS) values from heat-transfer-and-solidification-only simulations of single-track laser scan on bulk IN625 to identify that an appropriate range for the model constant “A” in Hunt’s model was between 1.3 and 1.7. Using the following values for the material properties of IN625 – k = 0.48, \( \Gamma \) = 2.2e-7 K m; \( \Delta T_{\text{o}} \)=60 K, and DL = 3.0e-9 m2/s—the factor \( \left( {k\Gamma \Delta T_{\text{o}} D_{\text{L}} } \right)^{0.25} \), was estimated to be 3.7131e-4 \( \left( {{\text{K}}^{2} {\text{m}}^{3} {\text{s}}^{ - 1} } \right)^{0.25} \). Based on Ghosh et al.,[35] considering the mean value of “A” of 1.5 for Hunt’s model constant, the constant \( A_{\text{o}} \) in Equation 11 was evaluated and PDAS equation becomes \( \lambda_{1} \left( {G, V} \right)\left[ {\mu {\text{m}}} \right] = 557\;G\left[ {\text{K/m}} \right]^{ - 0.5} V\left[ {\text{m/s}} \right]^{ - 0.25} \). Based on the review presented in this Section, it was decided not to include Kurz–Fisher model in this study. The PDAS formulations included in this study are shown in Table I, as: (a) original Hunt’s model (i.e., with A = 2.83 and \( A_{\text{o}} \) = 1051), (b) adjusted Hunt’s model by Ghosh et al.,[35] in which HTS-only data was used in the model parameter identification, (c) adjusted Hunt’s model using HTS-FD data in this study, and (d) Trivedi–Kurz (TK) model.

Table I PDAS Correlations Used in This Study
Full size table

Experimental Data

For IN625, the melt-pool shape was found to be similar for both powder and fully-dense substrate cases, although the LPBF exhibited additional width and height variations due to handling of the powder particles.[35] Thus, experiments for single-track laser fusion on fully-dense bulk material are a precursor to the LPBF. Moreover, the microstructure for single-track laser is representative for the LPBF process as the bottom part of the melt-pool would consist of re-melted and re-solidified bulk material, which is below the powder-bed, resulting in similar solidification conditions as those for STLF. The laser beam diameter, db, was 80 μm.

In order to assess the fluid dynamics effects on the microstructure, a series of four STLF cases were selected (Table II) with different primary dendrite arm spacings (PDAS). The high-resolution SEM micrographs for the melt-pool are shown in Figure 1. These SEM pictures were sized to scale, i.e., the 20 micron marker has the same length in all the four figures. The laser processing conditions are indicated on the top left corner of teach SEM picture; the laser power (W) and speed (mm/s) are indicated following the letter “p” and “s,” respectively. The extent of the meltpool is marked with a solid curved line while the extent of the equiaxed region in the center of the meltpool is marked with a dotted line. On the vertical centerline of the meltpool, the depth of the equiaxed domain was approximately 25, 10, 10, and 5 μm for p1s1, p1s2, p1s3, and p2s4, respectively.

Table II Cases Considered for Melt-Pool Assessment and Microstructure Analysis
Full size table
Fig. 1
figure 1

SEM micrographs of the melt-pool in bulk material for single laser track experiments: (a) p1s1, (b) p1s2, (c) p1s3, and (d) p2s4

Full size image

In Figure 2, the PDAS measurement locations are identified with thin and solid lines in a close-up view of the SEM micrograph for each of the cases considered. The PDAS was measured based on the linear intercept method, as: \( \lambda_{1} = \frac{L}{n - 1} \), where \( \lambda_{1} \) is the PADS, \( L \) is the length of the sampling line, and \( n \) is the number of dendrite arms that intersects the sampling line. In Figure 2, each sampling line is shown with a thin-solid line and each corresponding PDAS value is shown in \( \mu {\text{m}} \) at its measurement location. The height of the domain where PDAS was measured from the bottom of the meltpool by is approximately, 38, 39, 25, and 19 μm for the p1s1, p1s2, p1s3, and p2s4, respectively (i.e., by excluding the top region of 20, 10, 15, and 12 μm, respectively). Using the individual PDAS measurements shown in Figure 2, the mean PDAS values and its standard deviation, σPDAS, for each case were calculated and are shown in Table II.

Fig. 2
figure 2

SEM micrographs showing locations where PDAS was measured: (a) p1s1, (b) p1s2, (c) p1s3, and (d) p2s4

Full size image

Setup of STLF Simulation Model and Material Properties

In this study, a single laser scan line was simulated until uniform steady-state conditions were reached. Two types of models are considered for the simulations: heat-transfer-solidification-only (HTS), and a coupled heat-transfer-solidification and fluid dynamics (HTS_FD), as illustrated in Table III. The interdendritic flow and recoil pressure were not included in the HTS_FD model. Using symmetry on the centerline of the scanline, the computational domain was set to 1.1, 0.5, and 0.4 mm in the scan direction (X), lateral direction (Y), and depth (Z), respectively. The mesh used is shown in Appendix A. Within the overall domain, a subdomain in which the melt-pool is expected to develop was further refined. The length of the high-resolution subdomain was of 0.709, 0.103, and 0.072 mm in the X, Y, and Z directions (Appendix A). The mesh resolution in the high-resolution subdomain was 4.07, 2.65, and 1.85 μm in X, Y, Z directions, respectively. The heat transfer boundary conditions for the computational domain are indicated in Figure 3.

Table III Simulation Types to Assess Fluid Dynamics Effects on STLF Simulations
Full size table
Fig. 3
figure 3

Computational domain for LPBFAM simulations identifying heat transfer boundary conditions

Full size image

Thermophysical properties for the numerical simulations were selected based on the literature review or were calculated based on thermodynamic simulations using JMatPro software.[52] The latent heat of solidification was 290 kJ/kg. The fraction solid was considered to vary linearly with temperature within the solidification range of [1290:1350] °C. Thermophysical property data is given Appendix B for specific heat, thermal conductivity, viscosity, and surface tension. There was no experimental data found in the literature on coefficient of surface tension specifically for IN625. The JMatPro calculated data on coefficient of surface tension, dσ/dT, was found to vary between − 4e-4 N/m at liquidus to − 2.4e-4 N/m at temperatures above 2500 °C (Appendix B). The experimental data for dσ/dT of Ni superalloys (CMSX-4, IN738LC, MM247LC, and C263) was found to vary between − 6.8e-4 and − 15e-4 N/m.[53,54] For IN718, dσ/dT was measured to be − 1.1e-4 N/m.[55] This range of reported values for dσ/dT are fairly close to those estimated based on JMatPro simulation data.

Two optical properties are required for the simulation: total hemispherical emissivity for thermal radiation losses and absorptivity at the laser wavelength, \( A_{\lambda } \). In this study, an emissivity of 0.7 was considered for thermal radiation losses. Keller et al.[51] used \( A_{\lambda } \) = 0.5 for modeling the LPBF of IN625 alloy. For an Inconel alloy (79.5 Ni, 13.0 Cr, 6.5 Fe, and 0.08 C), the absorptivity at 1.06 μm wavelength was estimated using a validated Drude model[56] to be approximately 0.24 at temperatures of 120 °C to 600 °C.[57] Using the following relationship between the absorptivity and resistivity r(T),[58,59]:

$$ A_{\lambda } \left( T \right) = 0.365\sqrt {\frac{r\left( T \right)}{\lambda }} - 0.0667\frac{r\left( T \right)}{\lambda } + 0.006\left( {\frac{r\left( T \right)}{\lambda }} \right)^{1.5} , $$
(12)

the absorptivity for IN625 was estimated to be 0.33 at 1.064 μm wavelength and an electrical resistivity of 1.3e-4 Ohm-cm.

Numerical Simulation Results for Melt-Pool Shape

To assess the fluid dynamics effects on the melt-pool geometry, numerical simulations were carried out for one case with a laser power of 100 W and scan speed of 300 mm/s. For this case, the experimental width and depth of the melt-pool were 132 and 46 μm, respectively.

First, the melt-pool dependence on evaporation and surface tension coefficient was studied for a laser absorptivity of Aλ = 0.5, to enable a direct comparison with results from other studies such as Reference 48. These STLF simulations were conducted on OLCF Titan using 288 CPUs per each run. Each STLF simulation ran for 12 hours to simulate approximately 2 ms of actual STLF processing time. Second, the melt-pool dependence was studied as a function of evaporation and surface tension coefficient for an absorptivity of Aλ = 0.3, value which, to the best of our knowledge, is believed to be more accurate for IN625.

Sensitivity of melt-pool geometry on evaporation, surface tension coefficient, and laser absorptivity

Five cases were considered for HTS_FD simulations for different values of the surface tension coefficient, dσ/dT, and laser absorptivity, Aλ (Table IV). The calculated melt-pool width, Wcalc, and melt-pool depth, Hcalc, are shown in Table IV. The corresponding errors of the calculated melt-pool widths and depths with respect to the measured values (Wexp and Hexp) were calculated, as:

Table IV Sensitivity of Melt-Pool Geometry on the Surface Tension Coefficient, dσ/dT, and Laser Absorptivity Using HTS_FD Model with β = 1 (P = 100 W; U = 300 mm/s)
Full size table
$$ W_{\text{error}} \left[ {\text{pct}} \right] = 100*\left( {W_{\exp } - W_{\text{calc}} } \right)/W_{ \exp } ,\;{\text{and}}\;H_{\text{error}} \left[ {\text{pct}} \right] = 100*\left( {H_{\exp } - H_{\text{calc}} } \right)/H_{\exp } $$
(13)

The errors for the two melt-pool dimensions were shown in the last two columns of Table IV. By comparing the melt-pool data at constant dσ/dT, an increase in absorptivity was found to slightly increase the melt-pool depth and width. The error in the melt-pool depth is high as compared to that for the melt-pool width. The error between computed and experimental values for the melt-pool width is lowest for an absorptivity values of 0.45 and 0.5. For the melt-pool depth, the lowest for an absorptivity value 0.55. Note that the results shown in Table IV were obtained with the maximum evaporation rate, i.e., for β = 1 (the empirical coefficient that accounts for condensation effects). An overall error between the calculated and measured melt-pool dimensions was simply defined as the average of the absolute values of the errors for the pool width and pool height, i.e., (abs(Herr) + abs(Werr))/2. The overall error for cases shown in Table IV was found to range from 17.9 to 27.6 pct, with the smallest overall error for the dσ/dT = − 2e-4 N/m and Aλ = 0.5.

In an attempt to rule out the fact that the discrepancy might be due to an inappropriate evaporation flux, the empirical coefficient that accounts for condensation effects, β, was varied as shown in Table V from its maximum value of 1 (for the “a” cases), to 0.2, 0.1, and 0.05. The pool depth was found to exhibit a weak dependence on the evaporation factor, β. The overall error, (abs(Herr) + abs(Werr))/2, was found to range from 20 to 29.4 pct, with the smallest overall error for β = 1. Note that for the simulation of laser processing and welding of superalloys a factor β = 0.1 was found to be appropriate.[60,61,62]

Table V Sensitivity of Melt-Pool Geometry on Evaporation Flux (Variable β) for dσ/dT = − 3e-4 N/m and Aλ = 0.5 Using HTS_FD Model (P = 100 W; U = 300 mm/s)
Full size table

Sensitivity of melt-pool geometry on evaporation and surface tension coefficient for a laser absorptivity of 0.3

For each of the p1s2, p1s3, and p2s4 cases, four HTS runs and three HTS_FD runs were conducted. Since the case p1s1 was for a very low laser scan speed, which required a much larger simulation time, only one HTS_FD case and three HTS cases were run. The calculated melt-pool width (in μm) and melt-pool depth (in μm) are shown in Tables VI, VII, and VIII. The corresponding errors between the calculated melt-pool widths and depths were shown in the last columns of these tables.

Table VI Calculated Melt-Pool Width and Height for Two Evaporation Fluxes (β = 0.1 and 0.5) Using the HTS Model (dσ/dT = − 3e-4 N/m) and Aλ = 0.5
Full size table
Table VII Calculated Melt-Pool Width and Height Using the Fluid Flow HTS_FD Model (β = 0.1, dσ/dT = − 3e-4 N/m) and Aλ = 0.5
Full size table
Table VIII Calculated Melt-Pool width and Height and β = 0.1 Using the Fluid Flow HTS_FD Model (dσ/dT = − 0.75e-4 N/m) and Aλ = 0.5
Full size table

Two values of 0.5 and 0.1 were selected for the empirical coefficient that accounts for condensation effects, β (Table VI). For cases at low laser power of 75 W, both the pool depth and pool width change with the evaporation flux. With negative surface tension coefficient, dσ/dT, it is expected that the calculated values from fluid flow simulations would result in increased pool widths and decreased melt-pool depths. For β = 0.5, the calculated pool depths were smaller than that the experimentally measured values, and the calculated values are expected to further decrease for the fluid flow simulations. Since β = 0.5 is unlikely to yield more accurate results for the fluid flow cases, cases with β = 0.5 were further excluded from further evaluation for the HTS_FD model.

The temperature coefficient of the surface tension, dσ/dT, was then decreased from − 3e-4, to 1.5e-4, and − 0.75e-4. Results for dσ/dT = − 0.75e-4 N/mK are shown in Table VIII and compared to those in Table VII (dσ/dT = − 0.75e-4 N/mK) . The decrease in the dσ/dT, from − 3e-4 to − 0.75e-4, was found to lower the error in the pool depth to levels below 10 pct while the pool width error further decreased by approximately 6 to 12 pct (Tables VII and VIII). An overall error between the calculated and measured melt-pool dimensions can be simply defined as the average of the absolute values of the errors for the pool width and pool height. The overall error was found to range from 23.6 to 32.1 pct for the results in Table VII (dσ/dT = − 3e-4 N/m) and from 12.7 to 17.4 pct for the results in Table VIII (dσ/dT = − 0.75e-4 N/m).

Since the microstructural unit of the LPBF parts is a single laser melted track, it is worth the study of a single-track, keeping in mind the re-melting of the top bulk layer of material and the re-melting of one corner of the melt-pool due to the overlap during the laser scan for the adjacent layer (Figure 4). Thus, the analysis of the microstructure distribution within a single-track is a focus of this study.

Fig. 4
figure 4

Schematic showing the representative microstructure region within a vertical cross-section normal to the scan direction of a single-track melt-pool. PDAS measurement region was also indicated

Full size image

Numerical Simulation Results for Microstructure-Related Variables

From numerical simulation results, thermal gradient, G, and solidification velocity, V, were obtained to predict the microstructure type. The spatial distribution of the microstructure variables in a vertical cross-section through the solidified track of liquid pool, normal to the laser scanning direction, were presented for p1s2, p1s3, and p2s4 cases. These simulations were run with the parameters identified in the previous section, as: absorptivity Aλ = 0.3, empirical coefficient for condensation effects β = 0.1, and temperature coefficient of the surface tension, dσ/dT = − 0.75e-4 N/m.

Distribution of microstructure-related variables within the melt-pool

For these three cases, the distributions of the thermal gradient, G, solidification velocity, V, and cooling rate, GV, are shown in Figures 5, 6, and 7, respectively. The same length-scale was used in the figures. The figure legends indicate the minimum and maximum values for each case to provide the best visual illustration of each variable distribution. As shown in Figure 5, the thermal gradient, G, was found to exhibit a maximum around the edge of the solidified pool and a minimum in the entire central region of the solidified pool for the heat-transfer-only (HTS) model (Figures 5(a), (c), and (e)). For the fluid dynamics (HTS_FD) model (Figures 5(b), (d), and (f)), G was found to exhibit a maximum at the extremity of the solidified pool (at the free surface and not over the entire edge of the solidified pool) and the minimum in the central region below the free surface of the solidified pool (and not at the free surface of the central region).

Fig. 5
figure 5

Distribution of thermal gradient, G, for the: (a) p1s2 (HTS), (b) p1s2 (HTS_FD), (c) p1s3 (HTS), (d) p1s3 (HTS_FD), (e) p2s4 (HTS), (f) p2s4 (HTS_FD)

Full size image
Fig. 6
figure 6

Distribution of solidification velocity, V, for the: (a) p1s2 (HTS), (b) p1s2 (HTS_FD), (c) p1s3 (HTS), (d) p1s3 (HTS_FD), (e) p2s4 (HTS), (f) p2s4 (HTS_FD)

Full size image
Fig. 7
figure 7

Distribution of cooling rate, GV, for the: (a) p1s2 (HTS), (b) p1s2 (HTS_FD), (c) p1s3 (HTS), (d) p1s3 (HTS_FD), (e) p2s4 (HTS), (f) p2s4 (HTS_FD)

Full size image

As shown in Figure 6 for the two cases at 75 W laser power, the solidification velocity, V, was found to exhibit a minimum around the edge of the solidified pool and a maximum in the central region of the solidified pool for the heat-transfer-only (HTS) model (Figures 6(a) and (b)). For the fluid dynamics model (HTS_FD), V was found to exhibit some local variations within the melt-pool. At a higher power (case p2s4), V was found to exhibit similar distributions for both HTS and HTS_FD models (Figures 6(e) and (f)).

As shown in Figure 7, the cooling rate was found to exhibit different distributions for the HTS model for each of the three cases considered (Figures 7(a), (c), and (e)), exhibiting a maximum in the melt-pool center for p1s2 case, two-peak value regions (melt-pool center and edge) for p1s3 case, and at the top surface edge of the melt-pool for p2s4 case. For the HTS_FD model at power of 75 W (i.e., cases p1s2 and p1s3), the maxima were located at different locations than those for the HTS model runs. For the p2s4 case, the distribution of the cooling rate is very similar for the HTS and HTS_FD models (Figures 7(e) and (f)).

The minimum and maximum values for G, V, and GV, over the entire melt-pool, were summarized in Table IX. The minimum G values for the HTS and HTS_FD simulations are fairly close, with slightly lower values attained for the fluid dynamic cases (HTS_FD). On the other hand, the maximum G values were observed for the HTS_FD simulations (with the exception of p2s4 case). Both the minimum and maximum values for V were found to occur for the HTS_FD simulations, with the exception of the p2s4 case in which the values were very similar between the HTS and HTS_FD simulations. For the HTS_FD simulations, the minimum values of the cooling rate, GV, were found to be approximately half than their corresponding values for HTS simulations. By contrast with the min GV values, the maximum values for GV were found to occur for the HTS_FD simulations. Thus, HTS_FD simulations were found to exhibit a wider range of cooling rates than the HTS simulations.

Table IX Minimum and Maximum Values for G, V, and GV Within the Melt-Pool for HTS and HTS_FD Models
Full size table

Analysis of Microstructure Distribution

Figure 8 indicates the location of vertical lines and horizontal lines along which the calculated data for microstructure-related parameters (G and V) was obtained and analyzed in this section.

Fig. 8
figure 8

Schematic showing the location of lines in the vertical and horizontal directions along which microstructure-related variables were obtained in a vertical cross-section normal to the scan direction

Full size image

The solidification map data shown in Figure 9, which was obtained along one vertical line, c, and four horizontal lines (t, h1, h2, and h3), for the HTS and HTS_FD models. The solidification map data indicates the following: (a) the highest solidification velocity and lowest thermal gradient is located near the top surface region in the center of the melt-pool, (b) the smallest solidification velocity is located always near the bottom of the melt-pool in the central region (“c” dataline). The results for the HTS_FD model shows less variation of the solidification velocity with the thermal gradient (i.e., flatter profiles) and exhibited an increased spread in the V(G) variation within the melt-pool with respect to the HTS-only model results (see Figures 9(a) and (b); (c) and (d); (e) and (f), respectively). This can be explained by the fact that the fluid flow convection has an effect of reducing the temperature variation within the meltpool. Among the three cases considered, the p2s4 case exhibits similar results for both HTS and HTS_FD models. This indicates that the influence from the fluid flow is not significant and has less impact on microstructure at low linear power density of the laser (i.e., Power/Speed).

Fig. 9
figure 9

Solidification maps V(G) for: (a) p1s2 (HTS), (b) p1s2 (HTS_FD), (c) p1s3 (HTS), (d) p1s3 (HTS_FD), (e) p2s4 (HTS), (f) p2s4 (HTS_FD)

Full size image

For p1s2, p1s3, and p2s4, Figures 10 and 11 show the spatial PDAS variation along the centerline (c-dataline) for Trivedi-Kurz (TK) model and H_HTS_fit model,[35] respectively. For PDAS for the HTS and HTS_FD models, were obtained as follows λΗΤS(GHTS, VHTS) and λFD(GHTS_FD, VHTS_FD), respectively. For the p1s2 case (Figure 10(a)), the PDASFD shows smaller values near the top surface, consistent with the smaller PDAS evidenced in SEM micrograph (Figure 2(b)), while the PDASHTS shows a monotonic decrease, with a maximum at the top surface. For p1s3, both PDASFD and PDASHTS exhibit an unrealistic behavior as compared to the corresponding SEM micrograph (Figure 2(c)), i.e., with the maximum values at the top surface. For p2s4, both PDASFD and PDASHTS exhibit similar behavior as G and V did not exhibit significant variation with the fluid flow. Unlike the other two cases discussed, the equiaxed region only extends for less than 5 μm (Figure 2(d)), and PDAS at the top surface is not as refined as that for p1s2 and p1s3 cases. Thus, the actual variation with the maximum PDAS at the very top surface is not considered to be unrealistic. With current model parameters, TK model clearly overpredicts the PDAS, as the mean measured PDAS values (Table II) were 0.7, 0.59, and 0.45 μm, for the p1s2, p1s3, and p2s4, respectively.

Fig. 10
figure 10

Variation of PDAS along the vertical centerline for TK model (Appendix C—Trivedi and Kurz[49]), for: (a) p1s2, (b) p1s3, and (c) p2s4

Full size image
Fig. 11
figure 11

Variation of PDAS along the vertical centerline for H_HTS_fit, \( \lambda_{1} \left[ {\mu {\text{m}}} \right] = 557\;G\left[ {\text{K/m}} \right]^{ - 0.5} V\left[ {\text{m/s}} \right]^{ - 0.25} \), based on Ghosh et al.[35] for: (a) p1s2, (b) p1s3, and (c) p2s4

Full size image

Figure 11 shows the spatial PDAS distribution for the H_HTS formulation, obtained with the adjusted Hunt’s model[35] as indicated in Table I. The PDAS data with the H_HTS model indicate the following: (a) for low laser powers (cases p1s2 and p1s3), the PDAS obtained with the fluid dynamics model (HTS_FD) is larger than that with the simple HTS model, (b) for p1s2, PDAS with the fluid dynamics model is larger by more than 35 pct as compared with the PDAS obtained with the HTS model (Figure 11(a)), (c) for p1s3, PDAS with the fluid dynamics model is larger by more than 30 pct as compared with the PDAS obtained with the HTS model (Figure 11(b)). For p2s4, PDAS exhibit similar values with and without the fluid flow as G and V did not exhibit significant variation with the fluid flow.

To quantify the quality of the agreement between the measured values for PDAS, λ1,exp, and the computed values, λ1,comp, the mean PDAS, λ1, were used to calculate the deviation of the mean PDAS, Δλ, for each case, as follows:

$$ \Delta \lambda_{{\text{HTS},\;m}} \left[ {\mu m} \right] = \left( {\lambda_{1,\;comp,\;m} \left( {G_{\text{HTS}} ,V_{\text{HTS}} } \right) - \lambda_{1,\exp ,\;m} } \right),\;{\text{and}}\;\Delta \lambda_{{\text{HTS},\;m}} \left[ {\text{pct}} \right] = 100\Delta \lambda_{{{\text{HTS}},\;m}} /\lambda_{1,\exp ,\;m} , $$
(14)
$$ \Delta \lambda_{{{\text{FD,}}\;m}} \left[ {\mu m} \right] = \left( {\lambda_{{ 1 , {\text{comp,}}\;m}} \left( {G_{{{\text{HTS\_FD}}}} ,\;V_{{{\text{HTS\_FD}}}} } \right) - \lambda_{1,\exp ,\;m} } \right),\;{\text{and}}\;\Delta \lambda_{{{\text{FD,}}\;m}} \left[ {\text{pct}} \right] = 100\Delta \lambda_{{{\text{FD,}}\;m}} /\lambda_{1,\exp ,\;m} , $$
(15)

where index m specifies the case id (p1s2, p1s3, or p2s4) and the average λ1,comp, m(GHTS, VHTS) and λ1,comp, m(GHTS_FD, VHTS_FD) were calculated using the data in the corresponding region where experimental data shown in Table I was measured for each case.

Using the PDAS deviations obtained for each case, the overall PDAS deviations ΔλΗΤS and ΔλFD were calculated as a simple average among the three cases considered for the models considered. The deviation for the HTS and HTS_FD models in order to assess the effect of the fluid dynamics and are shown in Table X. The deviation of the PDASTK was found to be approximately 110 pct (Table X). An attempt was made to enhance the TK model by adjusting the only model parameter that is not related to a material property. The data from the three cases obtained with the HTS_FD TK model, i.e., λp1s2(GHTS_FD, VHTS_FD), λp1s3(GHTS_FD, VHTS_FD), and λp1s4(GHTS_FD, VHTS_FD),was curve fit with a least-square procedure to minimize the corresponding measured values λexp-p1s2, λexp-p1s3, and λexp-p2s4, respectively. This effectively consists obtaining the actual multiplier of the PDAS data shown in Figure 10 that would yield the best fit with experimental data. This adjusted TK model is denoted TK_FD_fit. The best fit with the experimental data yields a deviation of approximately 19 pct (Table X). The deviations of the PDAS values for the H_HTS PDAS model, which is the adjusted Hunt’s model,[35] were found to be approximately − 40 and – 30 pct for the HTS and HTS_FD process model (Table X). Simply running the original Hunt’s PDAS model yields overpredictions in the mean PDAS of approximately 11 and 27 pct, for the HTS and HTS_FD process model. Similar to the TK_FD_fit for the TK model, the HTS_FD results for (G, V) were used to adjust the Hunt’s model. The deviations of the PDAS for the FD PDAS model, were found to be the lowest among all the PDAS models considered at approximately 1 pct (Table X). To further assess the applicability of the FD model, the average PDAS was calculated using data obtained along the three vertical datalines indicated in Figure 8 and the along horizontal dataline (h3). For a direct comparison with the measured PDAS, the dataset was confined to the central region just above the meltpool, which is similar to the region where experimental data were measured (Figure 4). This average PDASFD is shown with the standard deviation for all the three cases considered in Figure 12. The data shows an excellent agreement between the experimental and calculated values for the cases considered.

Table X Deviation Between Average Values Computed Using the HTS and HTS_FD Models with Respect to the Measured Average PDAS Values
Full size table
Fig. 12
figure 12

Mean PDAS for all three cases considered using the FD model \( \lambda_{1} \left[ {\mu {\text{m}}} \right] = 832\;G\left[ {\text{K/m}} \right]^{ - 0.5} V\left[ {\text{m/s}} \right]^{ - 0.25} \)

Full size image

Conclusions

In this study, multi-physics simulations for single-track laser fusion of IN625 alloy were conducted using a highly parallel open-source code, Truchas, considering a flat-free surface of the molten alloy, heat transfer, phase-change, evaporation, and surface tension phenomena. Results from numerical simulations were compared with experimental data on both meltpool shape and PDAS. Two sensitivity studies were conducted by varying the evaporation flux and the surface tension coefficient at fixed laser absorptivities of 0.5 and 0.3, respectively. The overall combined error between the calculated values with the HTS_FD model and measured values for both width and height of the melt-pool dimensions was attained for a surface tension coefficient of − 0.75e-4 N/m.

The spatial distribution of the solidification variables in a vertical cross-section through the solidified track of liquid pool, normal to the laser scanning direction, were obtained. For the fluid dynamics HTS_FD model, G was found to exhibit a maximum at the extremity of the solidified pool (i.e., at the free surface). By contrast, for HTS simulations, G was found to exhibit a maximum around the entire edge of the solidified pool. For the HTS_FD simulations, the minimum values of the cooling rate, GV, were found to be approximately half than their corresponding values for HTS simulations. By contrast with the min GV values, the maximum values for GV were found to occur for the HTS_FD simulations. Thus, HTS_FD simulations were found to exhibit a wider range of cooling rates than the HTS simulations, exhibited an increased spread in the V(G) variation and a decrease in the variation of solidification velocity with the thermal gradient (i.e., flatter profiles).within the melt-pool with respect to the HTS model results.

The solidification map data indicates that the highest solidification velocity and lowest thermal gradient is located near the top surface region in the center of the melt-pool. The minimum solidification velocity was found to be always located always near the bottom of the melt-pool in the central region. The results for the HTS_FD process model also show less variation in the solidification velocity with the thermal gradient (i.e., flatter profiles) and exhibited an increased spread in the V(G) variation within the melt-pool with respect to the HTS-only model results. This can be explained by the fact that the fluid flow convection has an effect of reducing the temperature variation within the meltpool.

Correlation models based on commonly applied power law dependence (i.e., Hunt’s model) and marginal stability theory (Trivedi–Kurz model) on thermal gradient and solidification velocity for primary dendrite arm spacing (PDAS) was extensively evaluated using the data on G and V from both HTS and HTS_FD simulations to quantify the effect of the fluid dynamics on the microstructure. The PDAS formulations included in this study are: (a) Trivedi-Kurz (TK) model, and (b) adjusted TK model using HTS-FD data in this study, (c) original Hunt’s model (i.e., with A = 2.83 and \( A_{\text{o}} \) = 1051), (d) adjusted Hunt’s model by Ghosh et al.,[35] in which HTS-only data was used in the model parameter identification, and (e) adjusted Hunt’s model using HTS-FD data in this study. The deviation of the PDASTK was found to be approximately 110 pct. The deviations of the PDAS values for the H_HTS PDAS model, which is the adjusted Hunt’s model, were found to be approximately − 40 and − 30 pct for the HTS and HTS_FD process model. Simply running the original Hunt’s PDAS model yields overpredictions in the mean PDAS of approximately 11 and 27 pct, for the HTS and HTS_FD process models, respectively. The deviations of the PDAS using the (G, V) results from the HTS_FD model, i.e., \( \lambda_{1} \left[ {\mu {\text{m}}} \right] = 832\;G\left[ {\text{K/m}} \right]^{ - 0.5} V\left[ {\text{m/s}} \right]^{ - 0.25} \), were found to be the lowest among all the PDAS models considered at approximately 1 pct.