Multi-Scale Modeling of Residual Stresses Evolution in Laser Powder Bed Fusion of Inconel 625

: Laser powder bed fusion exhibits many advantages for manufacturing complex geometries from hard to machine alloys such as IN625. However, a major drawback is the formation of high tensile residual stresses, and the complex relationship between the process parameters and the residual stresses has not been fully investigated. The current study presents multi-scale models to examine the variation of process parameters on melt pool dimensions, cyclic temperature evolutions, cooling rate, and cyclic stress generation and how they affect the stress end state. In addition, the effect of the same energy density, which is often overlooked, on the generated residual stresses is investigated. Multi-level validation is performed based on melt pool dimensions, temperature measurements with a two-color pyrometer, and ﬁnally, in-depth residual stress measurement. The results show that scan speed has the strongest effect on residual stresses, followed by laser power and hatch spacing. The results are explained in light of the non-linear temperature evolution, temperature gradient, and cooling rate during laser exposure, cooling time, and the rate during recoating time.


Introduction
Laser powder bed fusion (LPBF) [1] is an additive manufacturing technology that has gained momentum in recent years due to design freedom to manufacture highly complex and customizable geometries along with its capability of processing metals and ceramics [2].These advantages make LPBF desirable for many industries where it can be used, for example, to manufacture a conformally cooled mold for the tooling industry [3] or light-weighting or to reduce the number of assembled parts, which makes it desirable for the automotive and aerospace industries [4][5][6].Despite the many advantages of LPBF, it still has a few drawbacks that need to be overcome, such as porosity, dimensional accuracy, high surface roughness, and tensile residual stresses (RS) [7,8].The first three drawbacks have been extensively studied experimentally and controlled by optimizing the process parameters [9][10][11][12].In addition, dimensional accuracy and surface roughness were also rectified by post-machining [2,13,14].However, the surface and near-surface RS was mostly tensile, regardless of the process parameters employed [15][16][17][18], which require post-processing such as stress relieving or peening to negate these undesirable stresses.Tensile RS is formed in LPBF due to the high-temperature gradient generated due to the large heat input under the laser beam between the solidifying material and surrounding cool powder and substrate [19].Therefore, to reduce the undesirable tensile RS, the effects of process parameters on the thermal history and temperature gradients on its formation mechanism need to be understood.
Modeling offers insight into the thermal cycles, temperature gradients, and cooling rates (CR) during LPBF, affecting stress formation and microstructure.In addition, predicting melt pool dimensions can determine if sufficient overlap between laser tracks will happen, as it affects the resulting density.A vast array of work has been performed on the modeling of LPBF targeting different scales or modeling approaches, ranging from melt pool scale to part scale models.The first approach is modeling single tracks in LPBF, which can be modeled using finite element (FE) or computational fluid dynamics (CFD).For FE models, the studies aimed at predicting melt pool dimensions and temperature gradients along the melt pool width and length [20][21][22][23][24]. CFD models also simulate single tracks but with the added benefit of including fluid dynamics effects such as the Marangoni effect, which drives the convection inside the melt pool, affecting the melt pool shape and CR that controls the final microstructure [25][26][27][28].Both models were extended to simulate the effect of adjacent tracks on melt pool dimensions and temperature gradients; however, beyond that, there is little focus on the impact of process parameters [29][30][31].
Moreover, these models overlooked the effect of the temperature gradients, CR, and process parameters on stress formation during solidification.The only studies that have addressed this particular area were focused more on the effect of scan strategies on RS [32,33].The second approach used was part scale simulation of part deformation and showed an agreement between the predicted and the experimentally measured deformations or RS [34][35][36].
Despite the varying degrees of success of each of the previously mentioned modeling approaches, there are a few shortcomings or challenges regarding FE modeling of the LPBF process [37,38] that need to be addressed.The first issue is the proper representation of the powder layer and its apparent thermal conductivity.Several studies assume the powder thermal conductivity as a fraction of the bulk material properties using the porosity ratio of the powder layer [20,24,39,40].However, by examining a unit volume of the powder layer, it is found that it comprises solid powder particles and the voids between particles filled with inert gas.Hence, the heat transfer coefficient consists of solid thermal conductivity, inert stagnant gas thermal conductivity, and radiation between the powder particles.The second issue is applying a Gaussian surface heat flux to represent the laser beam heat input, adopted from welding modeling [32,[41][42][43].The problem with the surface heat flux is that it neglects the optical penetration of the incident laser through the voids between the powder particles and its consequent decay.
Conversely, a volume heat source that considers the powder layer's laser penetration relative to the powder size is essential to represent the heat input accurately.The third issue is the incapability of the FE models to include the Marangoni effects [43], which is important as convection is the dominant heat transfer mode inside the meltpool and will, in return, affect the melt pool dimensions and temperature.The fourth issue is the lack of experimentally measured laser absorptivity, which is important as the powder bed absorbs a portion of the incident laser beam energy.Several studies either assume the value for the laser absorptivity or use some analytical models to calculate a more accurate value [44,45].In addition, there is a lack of experimental temperature measurements for validation of the numerically predicted temperature [43,45].
The current work addresses the gaps in the literature regarding the modeling of LPBF of Inconel 625 using the following:

•
Multi-scale modeling, including single track, multi-track layer, and part scale models;

•
Comparing different volumetric heat sources with optical penetration;

•
A proper powder material model is used to calculate the temperature-dependent thermal conductivity incorporating the solid conductivity, the nitrogen gas conductivity, and the inter-particle radiation;

•
The Marangoni effects inside the melt pool by calculating the Marangoni number and its equivalent thermal conductivity;

•
Validation of the melt pool dimensions for single tracks with the literature; • Experimental temperature measurement using a two-color pyrometer to validate the temperature from the multi-track layer model;

•
In-depth RS profile measurement using XRD to validate the part scale model RS predictions; • Study the effect of process parameters on thermal history, temperature gradients, CR, and their subsequent impact on the formation of RS.

Transient Thermal Model
Heat transfer in powder bed fusion is governed by heat conduction defined by the transient heat transfer energy balance equation, Equation ( 1) where: k: Thermal conductivity of the powder bed (W/m•K) T: Temperature (K) .
q: Energy generated per unit volume (W/m 3 ) ρ: Density of the powder bed (kg/m 3 ) C: Specific heat capacity (J/kg•K) t: time (s) The initial condition for Equation ( 1) is set as the powder bed preheat temperature of 80 • C and represented by Equation ( 2)

Heat Flux Boundary Condition
As the laser beam irradiates the surface of the powder bed, heat is added into the powder, which can be represented as a specified heat flux boundary condition (BC) and defined using the Fourier heat flux equation given in Equation (3):

Convection Boundary Condition
A constant nitrogen gas flow is pumped into the build chamber to reduce the oxygen percentage and carry away any spatter and ejected particles from the melt pool.The gas flow over the powder bed surface draws away heat due to forced convection, which is calculated using Equation (4): where: h o : Overall heat transfer coefficient due to forced convection (W/m 2 •K) T pb : Powder bed surface temperature T ∞ : Ambient temperature (assumed 300 K) The heat transfer coefficient due to forced convection is calculated similarly to the flow over a flat plate.Therefore, h will depend on the speed of the inert gas, as in Equations ( 6) and (7), and all fluid flow properties are to be calculated at film temperature [32,46].The film temperature is the average temperature between the powder bed surface temperature and the gas flow temperature, Equation (5).The EOS M280 machine used to perform the experiments has an inert gas velocity of 3 m/s. where: T f : The film temperature (K) k g : Gas thermal conductivity (W/m Another source of heat loss is radiation heat transfer to the surroundings, which is defined using Equation (8): where: e: Emissivity of the powder bed Since the powder bed is porous, the bulk material properties cannot be used.Therefore, an analytical model is used to calculate the emissivity of the powder bed [47,48], represented by Equations ( 9)- (11). where: A H : Porous area fraction of powder surface ϕ: Powder bed porosity ratio e s : Emissivity of bulk material The powder bed porosity depends on the recoating system and powder flowability [9].However, it difficult to precisely measure the powder bed porosity and is usually assumed to be between 40% to 60%.

Material Model
In PBF, the powder is heated consistently up to melting.Therefore, the material model is governed mainly by the temperature-dependent thermophysical properties of the powder.An analytical model developed by Sih and Barlow [47] is used to calculate the thermal conductivity of the powder bed based on the bulk material thermal conductivity and powder bed porosity ratio.The model, Equation (12), considers the radiation between particles, Equation (13), and the presence of stagnant inert gas in the voids between the powder particles.
where: k pb : Powder bed thermal conductivity (W/m•K) k g : Inert gas thermal conductivity (temperature dependent) (W/m•K) ϕ: Powder bed porosity ratio k s : Thermal conductivity of the bulk material (temperature dependent) k r : Equivalent conductivity due to thermal radiation between powder particles e: Emissivity of the powder (calculated previously) T: Temperature (K) D p : Average diameter of the powder particles (m) The next two properties to model powder bed thermal behavior are the specific heat capacity and density.A simplistic linear is used based on the porosity of the powder bed The temperature-dependent thermo-physical properties of bulk IN625 were obtained from [49], and the full material properties used are shown in Table 1.The selection of a suitable material mode is essential to predict temperature history better.Since the primary mode of heat transfer in L-PBF is conduction, as mentioned above, a low thermal conductivity will lead to a high concentration of heat at the center of the melt pool and raise the predicted peak temperature.Employing the powder model shows that the solid thermal conductivity is 20 times the powder thermal conductivity, as shown in Figure 1.

Heat Source Modeling
A volumetric heat source model is employed to consider the laser beam penetration into the porous powder bed.Three models are chosen for comparison, the Goldak double ellipsoidal model, shown in Equation ( 16), which was developed for welding [52].
where P is the laser power (W), A c is the absorption coefficient, x, y, and z are the local, a, b, and c are the longitudinal, transverse, and depth dimensions of the heat source, respectively, v is the laser scan speed (mm/s), t is the time (s), and h is the hatch spacing (mm).The dimensions a and b were assumed as the radius of the laser spot, while the depth dimension c is set as the melt pool depth determined experimentally.
The second and third heat source models examined are categorized as attenuated volumetric heat sources, based on the assumption that, as the laser penetrates the inter-particles voids, it becomes scattered and loses its energy at a certain depth.This attenuation depth is called the optical penetration depth (OPD), and its value was measured for pure nickel [53] at different powder particle size distributions (PSD), as shown in Table 2.However, there are no measured values for OPD in IN625 powders, but since IN625 comprises 58% Ni, the values for pure nickel were used to calculate the OPD using interpolation.The experiments were performed using an IN625 powder with a PSD of − 45 µm [50], which gives an OPD of 100 µm.The difference between the two attenuated volumetric heat sources is the attenuation profile being either a linear decay model, represented by Equation (17), or an exponential decay model, calculated using Equation (18).
where P is the laser power (W), A c is the absorption coefficient, x, y, and z are the local, d is the OPD (mm), r is the laser spot radius (mm), v is the laser scan speed (mm/s), t is the time (s), and h is the hatch spacing (mm).

Melt Pool Modeling
The formation of a melt pool occurs in two steps: the first step occurs between the solidus and liquidus temperatures; hence, it is governed by the latent heat of fusion, given in Table 1.The second step is the fully melted phase, where fluid mechanics such as the Marangoni effect and natural convection dominate heat transfer in the melt pool.Marangoni effects are important, as they drive fluid flow inside the meltpool, hence altering its dimensions [54].Additionally, the Marangoni effect will influence the solidification process and the resulting grain formation [55,56].Simulation of the Marangoni effect requires the use of CFD; however, the current study employs an FE model to predict the RS.In order to overcome the absence of fluid flow modeling in FE, an effective thermal conductivity model [57] is used to account for heat transfer due to the Marangoni effect artificially, as shown in Equation (19).
where k e f f is the effective thermal conductivity (W/m•K), k l is the quiescent liquid thermal conductivity (W/m•K), h is the convective heat transfer coefficient (W/m 2 •K), and L is the characteristic length (m), which is assumed as half the melt pool width [57,58].The value of the convective heat transfer coefficient inside the melt pool can be calculated from the Nusselt number using Equation (20), which is the ratio of convective to conduction heat transfer.Since convection in the melt pool is driven by the Marangoni effect, then the Nusselt number can be calculated as a function of the Marangoni number shown in Equations ( 21) and (22).
where Nu is the Nusselt number, Ma is the dimensionless Marangoni number, dσ dT is the temperature gradient of surface tension (N•m/K), L is the characteristic length, ρ is the liquid density (kg/m 3 ), C pl is the liquid specific heat capacity, k l is the quiescent liquid thermal conductivity (W/m•K), µ is the dynamic viscosity (kg/m•s), and ∆T is the maximum temperature gradient in the melt pool, which is taken as the temperature difference between the center of the melt pool (peak temperature) to the boundary of the melt pool (solidus temperature) and is calculated using Equation ( 23)

Modeling Approach
The modeling approach is based on three models, each focusing on predicting a certain set of results to characterize the thermo-mechanical response occurring in L-PBF under different process parameters.The first model is a single-track model, with a very fine mesh, to test the validity of the proposed material, melt pool, and heat source models by predicting melt pool dimensions, temperature gradients, and CR and comparing them to experimental results.The second model is a layer scale model with a coarser mesh to simulate the effects of multi-tracks on the temperature history, CR, melt pool dimensions, and the subsequent RS evolution in a single layer.The third model is a low-fidelity part scale model, which examines the far-field temperature histories and how it drives the final in-depth RS and part distortion.The user subroutine USDFLD changes the material model from powder to liquid and then to solid as the material solidifies.The criteria for changing the material state depends on the condition of reaching the melting temperature, as shown below: The melt pool dimensions, temperature gradient, and CR are extracted from the model when the melt pool reaches a steady-state, determined when the peak temperature reaches a constant value.The process parameters studied include two laser powers, 140 W and 220 W, and two scan speeds, 500 mm/s and 650 mm/s.

Multi-Track Model
The multi-track model aims to observe the effect of hatch spacing, laser power, and scan speed on the temperature profile, temperature gradient, CR, and consequent effect on RS evolution during the solidification layer.In order to do so, a sequentially coupled thermal stress model is used, since, in LPBF, the temperature change will drive a change in stress, but the generation of stress will not cause an increase in temperature.Therefore, two separate models, 3D FE models, were built: firstly, a pure heat transfer model to calculate the temperature fields generated due to laser exposure.Secondly, the temperature fields from the heat transfer model were imported into a static stress model to calculate the thermal stresses developed due to the cooling down of the material during solidification.
The heat transfer model has a powder layer with dimensions of 4 mm × 4 mm × 0.04 mm and is tie-constrained to the substrate having a thickness ten times the powder layer.The model was meshed with DC3D8 hexahedral elements having a uniform size of 40 µm.In contrast, the substrate had a coarser mesh with the element size biased towards the top surface, where the element size starts at 40 µm and increases to 140 µm, as shown in Figure 3.The user subroutine USDFLD changed the powder layer's material properties from powder to liquid and solid material due to material heating, melting, and cooling.The scan strategy employed was a stripe scan with serpentine laser motion [50].The laser heat load was defined as a body flux acting on the powder layer, and the spatial motion of the laser beam was implemented using the DFLUX subroutine.The step time was set as the total time required to expose a layer with a cross-section of 10 mm × 10 mm, the same size as the experiments [50].The laser process parameters are stated in Table 3.A preheat temperature of 80 • C was assigned to the whole model and set as a BC for the bottom surface of the substrate.Convection and radiation are defined for the top surface of the powder layer with an ambient temperature of 40 • C. The temperature profiles, temperature gradients, and CRs are extracted at the center point of the powder layer to examine the effect of multi-track laser exposures and hatch spacing on temperature evolution.The static stress model had the same dimensions and mesh size as the heat transfer model.Thermal stresses are the main driver of RS and deformation in LPBF; hence, the temperature field from the heat transfer model was assigned as an initial condition in the stress model.Since LPBF is characterized by rapid cooling, the molten material solidifies once the laser beam passes.This means that solidification and stress formation progressively follow the laser beam during the multi-track exposure.Therefore, to model this behavior, "quiet" or "inactive" elements represented the powder material that did not develop any stresses yet.Once the laser beam passes over these elements, they are activated to calculate stress evolution during cooling.The user subroutine UEPACTIVATIONVOL was used for progressive activation of the powder layer elements in tandem with the laser scan strategy.The coordinates of the laser position are calculated at every increment, and the elements within a circle, with a diameter equal to the melt pool width calculated from the heat transfer model, are activated [59].Due to a limitation in Abaqus, a set of elements with an area equal to the laser spot diameter needed to be activated priory to avoid convergence issues with UEPACTIVATIONVOL.In order to reduce computational time and improve the convergence of the non-linear model, the unsymmetric stiffness matrix storage option is selected to use Newton's method to solve the non-linear problem.In addition, the strain state extrapolation feature is turned off to improve convergence [60,61].

Multi-Layer Model
A large-scale FE model is required to predict stresses and deformations on a part scale.However, utilizing a fine mesh from the high-fidelity model would be computationally infeasible.Therefore, a low fidelity model that predicts the far-field temperatures without focusing on the melt pool area during printing and recoating was proven to provide good predictions of part scale stresses and deformations [17,61].The ABAQUS AM module [62] was used to build the sequentially coupled temperature-stress part scale model, consisting of a heat transfer step followed by a stress analysis step, and finally removing the part from the build plate.
The heat transfer model consists of a 10 mm × 10 mm × 10 mm cube coupon meshed with DC3D8 200 µm elements, while the substrate has a coarser mesh of 2 mm in size, as shown in Figure 4.The 200 µm element used is a lump sum of five layers of 40 µm thickness to save on computational cost, since the interest is the far-field temperature history [17,61].The part and the substrate are connected using a "tie-constraint."Simulation of the PBF process requires successive addition of material using UEPACTIVATIONVOL subroutine based on an event series that specifies the time instances where the recoating starts and ends.A Matlab code was written to generate the recoater event series file, with the recoating time set at 5 s.The laser heat source is modeled as a point source, which is acceptable when the laser spot size is relatively smaller than the element size.A stripe scan with serpentine laser motion was used for in-layer exposure, while a 0 • -90 • , alternating scan was employed for subsequent layer exposures, as illustrated in Figure 5.An event series data set is created using Matlab to determine the time, spatial coordinates, laser power, and laser spot diameter corresponding to the employed scan strategy.The tool path-mesh intersection built-in module determines the intersection of the laser event series with the element integration points and then calls the user subroutine UMDFLUX to apply heat flux at the corresponding time [61,62].The newly formed external surfaces are automatically determined, and the radiation and convection heat transfer is applied.Finally, after the whole part is exposed, it is left to cool down for 600 s.The same geometry and mesh size were used for the stress analysis model, where the temperature fields from the heat transfer model are imported into the stress model.The thermal strains and stresses were calculated by imposing an initial temperature of 1000 • C, above which thermal stress is negligible.The difference between this temperature and the imported part temperature will determine the thermal stress field [17,61].Temperaturedependent material properties were used [49,63], while Johnson-Cook (J-C) plasticity [64] was used to model the temperature-dependent plasticity according to Equation ( 25): where σ is the von mises flow stress, ε is the equivalent plastic strain, .
ε is the equivalent plastic strain rate, and .ε o is the reference plastic strain rate.T is the current temperature, T m is the melting temperature, and T r is the reference temperature at which the J-C parameters were determined.A is the initial yield strength, B is the strain hardening coefficient, n is the strain hardening exponent, C is the strain rate coefficient, and m is the thermal softening exponent.The values for the J-C parameters used in the current study are listed in Table 4.The J-C parameters listed above were determined for a conventionally manufactured IN625; in an attempt to partially modify these parameters to match that of the additively manufactured IN625, a tensile test was performed at room temperature.A sub-size flat specimen was printed with dimensions according to the ASTM E8 standard [66] with a gauge length of 25 mm, a width of 6 mm, and a thickness of 5 mm.The true stress-strain curve was plotted, and the Ramberg-Osgood model [67] was used to fit the data according to Equation ( 26) and calculate the yield strength, the strength coefficient (K), and the strain hardening exponent (n), which are analogous to the A, B, and n parameters in the J-C equation.
A total fixation BC is applied to the bottom surface of the substrate.After the part was built and cooled down, the substrate was removed using the model change feature to simulate the part's separation from the substrate.The substrate removal will release some stresses, and the RS inside the part will remain.For validation, the stresses are extracted and averaged across an area 2 mm in diameter to resemble the spot size of the experimentally measured RS using XRD measurement.

Laser Absorptivity Measurement
Diffuse reflectance spectroscopy (DRS) was used to measure the optical absorptivity of the IN625 in a wavelength range of 400-1400 nm.The test measures the powder's reflectance, which is used to calculate the Kubelka-Munk absorption coefficient, as shown in [9].For example, for an EOSM280, the Yb-fiber laser used has a wavelength of 1070 nm that corresponds to an absorptivity coefficient of 0.62.

Single Track Validation
The melt pool width and depth predicted by the single-track model are validated by comparing them to the experimental melt pool dimensions measured in the literature [68,69].The single exposed by NIST [68] and Dilip et al. [69] were built with a layer thickness of 20 µm, and the model geometry was modified to accommodate the smaller layer thickness.The meltpool dimensions were measured by extracting the distance from the center of the melt pool to the outer temperature contour having a temperature corresponding the solidus temperature of 1290 • C.

Multi-Track Temperature Measurement
A Fluke Endurance-series two-color pyrometer was used to measure melt pool measurements during exposure of a single layer.The two-colour pyrometer operation is based on the ratio of two separate infrared bands with a slight difference in wavelength to determine the object's temperature.Hence, temperature measurement would be independent of the emissivity, thus reducing measurement uncertainty.The pyrometer used has a Si/Si detector with a nominal wavelength of 1 µm, and the ratio between the two wavelengths, also called the slope, was calibrated as 0.96 [70].The pyrometer has a temperature range of 1000-3000 • C, and it was installed inside the build chamber of an Omnitek L-PBF machine, as shown in Figure 6a, at a distance of 85 mm, which results in a measurement spot size of 2 mm in diameter.For validation, the elements temperatures above 1000 • C, i.e., the lower limit of the pyrometer, within a 2 mm diameter circle were averaged and compared to the pyrometer measurements [71], as shown in Figure 6b.The tested parameters were 220 W, 650 mm/s, 0.12 mm hatch spacing, and a layer thickness of 0.04 mm.Since the pyrometer had a response time of 10 ms, the temperatures calculated from the model are extracted and averaged at an increment of 10 ms.

Residual Stresses Measurements
In-depth RS were measured in IN625 cube coupons [50] using X-ray diffraction (XRD).Since X-ray penetration depth is about 5 µm in Inconel, electro-polishing was used to minimally affect the existing stress field [72].The material was removed in 0.1 mm increments down to a depth of 1 mm from the top surface.RS were measured in two directions; parallel to the laser motion (scan direction) and perpendicular to it (hatch direction).The in-depth RS measurements were performed on two samples corresponding to two laser powers, 140 W and 220 W, at a hatch spacing of 0.12 mm, and a constant scan speed of 650 mm/s.In addition, predicted RS is compared to surface RS measured by the current authors in a previous study [50] at different process parameters.

Single Tracks
The temperature contours of the predicted melt pool are shown in Figure 7 for the 20 µm layer thickness validation model with an exponential heat source.The grey-colored area represents temperatures above the solidus temperature of 1290 • C. The temperature profiles are plotted in the width direction, Figure 8, and in-depth, to calculate the melt pool dimensions, as shown in Figure 9.The intersection point of the temperature profiles with the solidus temperature locus determines the width and depth of the melt pool.Comparison of the experimentally measured melt pool width [68,69] with the predicted melt pool width for the three heat source models is shown in Figure 8.The exponential heat source model predicted a melt pool width of 131 µm, which best agrees with the experimental measurements.However, the melt pool depth was underestimated for all three heat source models, as shown in Figure 9.The melt pool depth measured by NIST [68] was much deeper than the predicted value due to printing on the bare substrate without a powder layer.Therefore, the exponential heat source was chosen for subsequent simulations.The length of the melt pool was measured in the study by NIST [68] using a high-speed infrared camera.By comparing the result with the numerical prediction, Figure 10 shows that the model greatly underestimates the melt pool length.The underestimation could be attributed to approximating the Marangoni effect inside the melt pool and calculating its value based only on the outward flow in the melt pool width direction.Calculating the melt pool width and depth for the parameters examined in the current study (Figure 11) shows that the melt pool dimensions decrease with increased scan speed, which agrees with the literature [68,69,73].Furthermore, it is observed that the melt pool width is inversely proportional to the square of the scan speed, while the melt pool depth is inversely proportional to the scan speed.
The temperature history is extracted at the midpoint along the laser scan path and plotted, as an example, for the case of 220 W and 650 mm/s, as shown in Figure 12.A rapid temperature rise is seen as the laser approaches the midpoint, followed by rapid cooling to 1290 • C. Finally, the CR is calculated as the slope of the temperature history profile once the temperature drops below the solidification temperature, which will drive RS.
The CR for all three scan speeds is plotted in Figure 13 and is to the order of 10 5 , characteristic of rapid cooling.The cooling increases linearly with increasing scan speed, regardless of the laser power.The temperature gradient is also plotted in Figure 13, as it affects the magnitude of thermal stresses formed during cooling according to Equation (27), where E is the modulus of elasticity (MPa), α is the thermal expansion coefficient, and ∆T is the temperature gradient (   The temperature gradient exhibits a downward trend with increased scan speed, resulting in RS with the same trend.However, previous experimental RS measurements show non-monotonic trends for RS values [50].These trends could mean that RS is governed by a temperature gradient, and the opposite effect of CR could play a role in determining the final stress state [42].Temperature predictions from the multi-track model are compared to pyrometer measurements, as shown in Figure 14.The numerical temperature predictions are in good agreement with the experimental data.However, the discrepancy at the highest temperature predictions could be due to underestimating the melt pool length; since the calculated temperature is averaged over the melt pool, a difference in the melt pool area will affect temperature prediction.The temperature history is plotted at three points, 1 mm apart, at 220 W laser power, 650 mm/s scan speed, and 0.12 mm hatch spacing, as shown in Figure 15.The three points exhibit the same thermal cycles where the temperature rises to a point shy of the melting temperature as the laser scans the precedent track.Afterward, the temperature at the point of interest cools down to about 500 • C as the laser beam passes directly over it, resulting in a peak temperature higher than the boiling temperature.Due to rapid cooling, the temperature drops from 2900 • C down to 570 • C, where the heat from the subsequent track scan causes the temperature to rise again to 950 • C. Beyond that point, the temperature drops to a plateau of about 650 • C until the whole layer is scanned.The CR and temperature gradients are plotted in Figure 16 for the points mentioned above.The presented CR is calculated as the temperature cools down below 1000 • C.Although a larger surrounding solid area could indicate that the cooling will increase because the solid material has higher thermal conductivity and specific heat capacity compared to the powder state, it is shown that the CR decreases as the area of the solidified layer increases.The reason is that the larger solidified volume retains more heat; therefore, the area surrounding the melt pool has a relatively high temperature, resulting in a lower temperature gradient, as shown in Figure 16, and lower heat transfer due to conduction.In addition, Figure 17 compares the temperature and von Mises stress evolutions during the layer exposure at different time instances until the end of exposure at 1.25 s.The peak temperatures are mainly focused inside and around the boundary of the melt pool and drastically decay beyond that.Conversely, the Mises stress is minimum around the melt pool boundaries and increase farther away from the melt pool due to the rapid cooling and stress generation, until it reaches its maximum value at the part edges, as seen in Figure 17f.The temperature evolution at the midpoint of the modeled area is plotted in Figure 18 for the total layer time of 6.25 s, corresponding to 1.25 s of laser exposure plus 5 s for powder recoating.The temperature profile can be divided into a rapid non-linear phase, which occurs during the first 1.25 s, followed by the linear phase.The thermal stresses are also plotted for the principal directions, the laser scan direction (S11), and the perpendicular direction, the hatch direction (S22), where S11 is observed to be always greater than S22.To understand why S11 is greater than S22, the RS formation mechanism needs to be explained.RS in LPBF is formed due to the difference in thermal strains during subsequent heating and cooling during subsequent scans called the temperature gradient mechanism (TGM) [15,73,74], governed by Equation ( 27).To better understand the formation mechanism of RS, Figure 19 shows a layer during laser exposure where the yellow area, A T , represents a track along the laser path surrounded by the previously solidified area, A S , and the subsequent powder area.This layer can be represented by the bar model analogy for simplification, whereas bar T and S are constrained at one end and free to move on the other end.As bar T is heated, it expands and pulls on bar S, resisting the motion by inducing compressive stress.The stress cycle undergone by the material is represented by the stress-temperature plot in Figure 19, where the compressive stress increases with temperature till reaching the yield temperature T y at point A, upon which plastic deformation occurs.As bar T cools down during solidification, the compressive strains and stresses decrease to zero at point C.As the temperature drops further, bar T shrinks and pulls on bar S, which resists this motion; therefore, bar T experiences tensile strains while bar S exhibits compressive strains.Finally, the tensile stress in bar T increases up to the yield point, where the material undergoes another plastic deformation until reaching equilibrium temperature.These plastic strains lock in the thermal stresses and prevent them from returning to zero, thus creating RS.The RS can be represented by Equation (28) or Equation ( 29), depending on the maximum temperature experienced by the material.
If T y ≤ T max < 2T y , then The same bar model analogy can be used to explain why S11 is greater than S22 by examining Figure 20, where the bar model is superimposed on arbitrary melt pool temperature contours [75].The thermal stress in the X direction (S11) is generated due to the temperature difference with the adjacent area in the Y direction (∆T Y ), while the thermal stress in Y direction (S22) arises due to the temperature difference in the X direction (∆T X ).As shown in Figure 16, (∆T Y ) is always greater than (∆T X ); hence, S11 is always greater than S22. Figure 21 provides a closer examination of the non-linear stress formation phase, where the material undergoes a few heating and cooling cycles.Commonly, stresses that are generated due to a moving heat source were considered to be affected by the temperature gradient only and its thermal softening effect on the material.However, in the material deforms at high strain rates due to the high CR, which causes strain rate hardening, thus opposing the thermal softening effect.Following the temperature and stress evolution in Figure 21, it is seen that the stress is zero as long as the temperature has not reached the melting temperature; hence, it remains in its powder state.As the laser passes over the point of interest during track n, the material melts, and as it solidifies, stresses are generated.By examining Figures 21 and 22 simultaneously, to observe the stress and plastic strain formation during cooling, it is seen that as the temperature decreases, the compressive plastic strain decreases, analogous to line BD in Figure 19, while tensile stresses increase.As the laser beam exposes track n + 1, the conducted heat causes a temperature rise of track n, as mentioned previously, thus causing an increase in compressive stress and a corresponding decrease in tensile stresses.Finally, the temperature due to the conducted heat reaches a maximum of 940 • C, the material cools down again, generating more tensile stresses.As the distance between track n and laser spot increases, the quantity of heat reaching track n decreases.Therefore, the material undergoes low peak heating and cooling cycles, which generates comparatively lower tensile stresses.

Effect of Scan Speed on RS
Examining the effect of the scan speed on the temperature history at the center point of the layer can be seen in Figure 23, where three scan speed levels of 500 mm/s, 650 mm/s, and 800 mm/s are shown.As expected, the temperature peaks decreased with increasing scan speed.Yet, the temperature at the layer midpoint did not reach the melting temperature during the exposure of track n − 1, as the hatch spacing was greater than the half-width of the melt pool.However, at a scan speed of 500 mm/s, the melt pool half-width was 120 µm, equal to the hatch spacing, but the peak temperature from track n − 1 was just shy of the melting temperature.Surface RS are measured after laser exposure plus a relaxation time equal to the recoating time of 5 s.It is shown in Figure 24 that tensile RS is formed, and its magnitude increases linearly with the increase in scan speed.The RS increase is due to the corresponding increase in the temperature gradient and CR, Figure 25, which would induce higher thermal stresses.The increase in the temperature gradient occurs even though increasing the scan speed will cause a reduction in the energy density added to the powder layer.The increasing trend can be attributed to the increase in CR, which causes the surrounding material to cool faster, thus creating a more significant temperature difference.In addition, as the speed increases, the melt pool width decreases, resulting in less heating of the adjacent tracks and lower temperature peak, as seen in Figure 23, where the amplitude of the third peak is lower at the higher scan speeds.Therefore, the lower temperatures of the adjacent tracks will lead to a larger temperature gradient.

Effect of Laser Power on RS
Three different laser power levels were modeled at a constant scan speed of 650 mm/s and 0.12 mm hatch spacing to examine the effect of laser power on the temperature profiles and RS.The temperature history in Figure 26 shows similar temperature evolutions for all three laser powers at the midpoint of track n.At any time, the minimum temperature is recorded with the 140 W power and increases with increasing the laser power.At 270 W, the melt pool width was 300 µm, which means that the half-width is larger than the hatch spacing.The larger inter-track overlap is reflected in the temperature history profile, where the point on track n reached the melting point while the laser was still scanning track n − 1, seen as the first peak on the red curve in Figure 26.Prediction of RS at the midpoint of track n shows that a laser power of 140 W resulted in the highest surface RS, as shown in Figure 27, and RS decreases with the increase in laser power, which agrees with the literature [50].The declining trend in RS is analogous to the temperature gradient and CR trend, shown in Figure 28, where the highest temperature gradient is recorded at 140 W.The decreasing trend can be traced back to the melt pool size, where the widths are measured to be 138 µm, 176 µm, and 298 µm, corresponding to 140 W, 220 W, and 270 W, respectively.Therefore, the wider melt pool will have a larger heat-affected zone (HAZ), as shown in Figure 29, and the surrounding material will have a higher temperature.The higher temperature of the HAZ is evident by the amplitude of the third peak in Figure 26, which increases with the increase in laser power.The higher the amplitude, the lower the temperature difference between the second (below solidus temperature) and third peak, corresponding to a lower temperature gradient.

Effect of Hatch Spacing on RS
The increase in hatch spacing from 0.08 mm to 0.12 mm resulted in a temperature history profile, Figure 30, different from that under varying scan speed or laser power.At 0.08 mm, the distance between subsequent tracks is small, so that the temperature of track n has peaked above the melting temperature three distinct times by the first three peaks on the blue curve in Figure 30.As the hatch spacing increases to 0.1, there are only two peaks above the melting point.At the largest hatch spacing, there is not enough heat conducted between tracks to raise the temperature of track n, thus leading to only one peak above the melting point.The increase in hatch spacing, i.e., the decrease in energy density, caused the surface RS to drop from 389 MPa to 307 MPa, as shown in Figure 31.Conversely, the decrease in energy density by increasing the scan speed or decreasing the laser power led to an increase in surface RS, as shown in Figures 24 and 27  It is found that the surface RS trend is governed by the temperature gradient, which had a more dominant effect on the surface RS than the CR, as shown in Figure 32.As the hatch spacing increases from 0.08 mm to 0.1 mm, the temperature gradient decreases; hence, the surface RS decrease.However, increasing the hatch spacing from 0.1 mm to 0.12 caused a reduction in the temperature gradient and an opposite increase in the CR, leading to a slightly lower surface RS.

Effect of Same Energy Density
Studying the effects of laser power, scan speed, and hatch spacing on surface RS shows opposing effects and, consequently, non-monotonic trends for the effect of the energy density on the surface RS.Therefore, the effect of the individual process parameters needs to be investigated at constant energy density.The first set of simulations are performed at an energy density of 3.67 J/mm 2 with a constant laser power of 220 W. Three levels for the scan speed, 500 mm/s, 600 mm/s, and 750 mm/s, and three hatch spacings, 0.12 mm, 0.1 mm, and 0.08 mm, are selected, respectively.
The temperature history, presented in Figure 33, shows that all three parameter combinations restulted in multiple temperatures peaks above the melting point.However, it is noticed that the maximum peak temperature was achieved at the lowest scan speed of 500 mm/s and the largest hatch spacing, which means that the scan speed has a more dominant effect on the maximum temperature versus the hatch spacing.Conversely, the hatch spacing has a more dominant effect when it comes to raising the temperature of the adjacent tracks, where the least hatch spacing resulted in two subsequent peaks with peak temperature highly above the melting point, as seen by the red curve in Figure 33.Examining the effect of the distinct temperature histories on the surface RS, presented in Figure 34, shows that the surface RS are lowest at 500 mm/s and 0.12 mm with a magnitude of 211 MPa and increase non-linearly to reach a maximum value of 454 MPa at 750 mm/s and 0.08 mm.The non-linear trend in the temperature gradient, shown in Figure 35, is the driver of the surface RS as the temperature gradient increases with the increase in scan speed and decrease in the hatch spacing.However, the CR decreases when the scan speed increases from 500 mm/s to 600 mm/s then increases beyond that point, which means that the temperature gradient mainly affects the surface RS.It can also be deduced that the scan speed has a more dominant effect than the hatch spacing on the surface RS, since increasing the hatch spacing from 0.1 mm to 0.12 resulted in a negligible change in the surface RS, as shown in Figure 31.However, the variation of the scan speed caused the surface RS to increase significantly.To further investigate the effect of process parameters on surface RS with the same energy density, the laser power and scan speed are varied, while the hatch spacing is kept constant at 0.12 mm.Three power levels, 170 W, 220 W, and 270 W, and three scan speeds, 500 mm/s, 650 mm/s, and 800 mm/s, gave an energy density of 2.82 J/mm 2 .All three parameter combinations resulted in almost identical temperature profiles, with a slight difference in the magnitude of the peaks, as shown in Figure 36.The results presented in Sections 4.2.2 and 4.2.3 show that the increase in scan speed caused an increase in surface RS; conversely, the increase in laser power caused a decrease in surface RS.However, the combined effect of increasing the laser power and scan speed increased surface RS, as shown in Figure 37. Therefore, the scan speed has a more dominant effect on the surface RS than the laser power.To investigate the underlying causes that drive this trend in the surface RS, the temperature gradients and CRs are plotted in Figure 38.The highest temperature gradient and CR were recorded at the lowest laser power and scan speed.Increasing the laser power and scan speed causes the temperature gradient and CR to decrease slightly and increase as the laser power and scan speed increase.Based on the previous observations, the surface RS followed the same trend as the temperature gradient.Therefore, if applied here, the surface RS should have been the highest at the 170 W/ 500 mm/s combination, decreasing at 220 W/650 mm/s, and then increasing again at 270 W/800 mm/s.However, this is not the case, as shown in Figure 37. Therefore, an in-depth investigation of the temperature and stress history is required.
The temperature and stress evolutions, plotted in Figure 39, show that, after the laser melts track n, the material starts to solidify and generates high thermal stresses.The stress reaches a magnitude of 460 MPa, 440 MPa, and 445 MPa, corresponding to 170 W/500 mm/s, 220 W/650 mm/s, and 270 W/800 mm/s, respectively.The trend of these stresses is similar to the temperature gradients shown in Figure 38.However, the final stress state is not reached yet, and further stress changes occur with the subsequent laser passes.As heat is conducted from the subsequent tracks, the temperature of track n rises, thus relieving the previously generated stresses causing the cyclic change in stresses shown in the left column in Figure 39.As the laser travels far enough from track n, the magnitude of the cyclic stress decreases.The conducted heat will keep reaching track n until the laser has finished scanning the whole layer; hence, the more time it takes for the laser to scan the layer, the more heat is conducted and, consequently, the greater the reduction in the stress magnitude.The laser scanning time from track n to the end of the layer is calculated at 0.84 s, 0.6 s, and 0.52 s, corresponding to 500 mm/s, 650 mm/s, and 800 mm/s, respectively, as shown in the right column of Figure 39.The difference in time durations allows the conducted heat to reduce the stresses to −80 MPa, −21 MPa, and 70 MPa, respectively.Once the laser finishes scanning, the temperature of the whole layer cools down for the same period, thus raising the stresses to their final state, where the lowest stress is recorded at 170 W and 500 mm/s.The J-C model's A, B, and n parameters were changed to the new value to compare the predicted RS in-depth profile with the experimental data.Due to the absence of tensile testing data for different process parameters, the same values of A, B, and n were used for all simulated process parameters.Additionally, since the measured Young's modulus is lower than its values reported in the literature, the newly measured value was used in the current study.Finally, Young's modulus was measured at room temperature.Still, Young's modulus data from the literature were offset by the difference between measured to obtain its values at high temperature.Their reported values were obtained at room temperature, as shown in Figure 41.

Numerical Results
The part scale model is used to predict the RS distribution across the whole part and the deflection of the part, which will affect its dimensional accuracy.The evolution of the Mises stress can be seen in Figure 42 at different time instances.At a time of 667 s, i.e., after a 4 mm height (Figure 42a), the maximum stresses are predicted at the bottom of the part, at the interface between the part and the substrate.As the number of layers increases to reach 10 mm, more heat is added to the top surface of the part and conducts downwards, causing the stress magnitude to decrease, as shown in Figure 42b.Once the laser exposure is completed, the part begins to cool down for 600 s, at which point the stresses in the part increase again (Figure 42c) to reach its maximum value.Once the part is separated from the substrate, the stress magnitudes decrease to reach their final value, as shown in Figure 42d.The RS were calculated at the end of the building process after a cooling time of 600 s and separation from the build plate.The RS were calculated by averaging the stresses across an area with a diameter of 2 mm.The predicted in-depth RS profile is validated by comparing to experimentally measured RS profile at two laser powers, 220 W and 140 W, a scan speed of 650 mm/s, and a hatch spacing of 0.12 mm, as shown in Figures 43 and 44   Two RS profiles are compared to the experimental results: the first RS profile is calculated using the J-C parameters found in the literature [65].The second RS profile is predicted with J-C parameters partially adjusted using data from tensile testing of IN625 manufactured using LPBF.It is found that the first profile under predicts the in-depth RS, while the second profile over predicts the RS profile.However, both profiles follow the same trend as the experimental results.These results highlight the need for J-C plasticity parameters for LPBF-manufactured IN625 to predict the generated stresses better.
Prediction of the in-depth RS with increasing the hatch spacing shows that the hatch spacing does not affect the RS profile, as shown in Figure 45, regardless of the laser power used.The hatch effect absence could be attributed to using a point source heat model, which does not generate the teardrop shape of a melt pool, so the effect of melt pool overlap is negated and would not affect the predicted RS.Examining the effect of laser power on the in-depth RS and deflection show that increasing the laser power from 140 W to 270 W decreases the in-depth RS, as seen in Figure 46.However, the generated RS for all three laser powers are highly tensile down to a depth of 1 mm from the surface.This trend is similar to that predicted from the layer model for the surface RS, as presented in Section 4.2.3.Increasing the scan speed also showed a similar trend to the results reported in Section 4.2.2.The in-depth RS increased with increasing the scan speed, as shown in Figure 47a.The deflection of the top surface in the Z direction is of particular interest, as it increases the chances of the recoater impacting the part and damaging or breaking the part off.Therefore, the deflection was extracted along the diagonal of the part, as shown in Figure 48a.It is observed that after cutting the part from the build plate, the edges of the part deflect upwards.At the same time, the center sinks in, resulting in a parabola-shaped surface.Examining the effect of process parameters on the deflection shows that the laser power has the same effect on the deflection as on RS, as seen in Figure 46b, where the lowest power resulted in the largest deflection.Conversely, the scan speed increased top surface deflection, as shown in Figure 47b.This is because the deflection of the part edges means that part of the generated stresses is released; hence, the magnitude of the stresses is always lower at the outer edges.

Conclusions
The current work presented a systematic approach to different modeling scales of laser powder bed fusion of IN625, including single tracks, multi-tracks, and multi-layer models.The shortcomings found in the literature were addressed by employing a Gaussian volumetric heat source model, a material model to accurately represent the thermosphysical properties of the powder layer, artificial inclusion of the Marangoni effect, and the use of the Johnson-Cook plasticity model to include the strain hardening, strain rate effect, and thermal softening of the material.The results from all models were validated with multiple experimental results, including melt pool dimensions, temperature measurements using a two-color pyrometer, and in-depth residual stress profiles.The main findings can be summarized as follows:

•
The exponential decay volumetric heat source, with an artificial increase in thermal conductivity to account for the Marangoni effects, resulted in the best prediction of melt pool dimensions.

•
The maximum residual stress is generated along the laser scan direction, governed by the temperature gradient and cooling rate in the hatch direction.

•
Increasing the scan speed resulted in an increase in the surface tensile residual stress due to the increase in the temperature gradient.

•
The surface tensile residual stress decreases with the increase in hatch spacing.

•
At the same energy density, the thermal stresses are mostly affected by the scan speed, laser power, and hatch spacing.

•
Although the high cooling rate will increase the strain rate experienced by the material, the evolution of the stresses is mainly dominated by the temperature gradient.

•
The cyclic temperature governs the final state of the stress during the non-linear phase and the cooling rate during the linear phase.

•
In-depth residual stresses on a part scale and the top surface deflection exhibit the same trend, where increasing the laser power lowers the residual stresses.

•
Increasing the scan speed increases the in-depth residual stresses as well as the top surface deflection.

2. 6 . 1 .
Single-Track ModelA 3D Lagrangian finite element (FE) model was built using ABAQUS to model singletrack exposure of IN625.The model dimensions are 1 mm × 240 µm × 300 µm, with the top 40 µm assigned as the powder layer, as shown in Figure2.The model was meshed with heat transfer DC3D8 elements with a minimum size of 10 µm employed in the powder layer.The volumetric heat flux was applied as a body on the powder layer using the user subroutine DFLUX to simulate the spatial movement of the laser beam with time.Symmetry BC was assigned to reduce the model's size, with the XZ plane taken as the plane of symmetry.

Figure 8 .
Figure 8. Validation of predicted melt pool width.

Figure 9 .
Figure 9. Validation of predicted melt pool depth.

Figure 10 .
Figure 10.Melt pool length prediction for different heat source models.

Figure 11 .
Figure 11.Single tracks melt pool dimensions at different process parameters.

Figure 12 .
Figure 12.Temperature history evolution for a single track.

Figure 13 .
Figure 13.Single tracks CR and temperature gradients.

Figure 17 .
Figure 17.Temperature contours evolution (a-c) and Mises stress evolution (d-f) at different time stamps.

Figure 18 .
Figure 18.Layer stress and temperature evolution at 220 W, 650 mm/s, and 0.12 mm.

Figure 19 .
Figure 19.RS formation mechanism with the bar model analogy.

Figure 23 .
Figure 23.Temperature history for different scan speeds at 220 W and 0.12 mm.

Figure 24 .
Figure 24.Effect of scan speed on the surface RS.

Figure 25 .
Figure 25.Effect of scan speed on (a) temperature gradient and (b) CR at 220 W and 0.12 mm.

Figure 26 .
Figure 26.Temperature history for different laser powers at 650 mm/s and 0.12 mm.

Figure 27 .
Figure 27.Effect of laser power on surface RS.

Figure 28 .
Figure 28.Effect of laser power on (a) temperature gradient and (b) CR at 650 mm/s and 0.12 mm.

Figure 30 .
Figure 30.Temperature history for different hatch spacing at 220 W and 650 mm/s. .

Figure 31 .
Figure 31.Effect of hatch spacing on surface RS.

Figure 32 .
Figure 32.Effect of hatch spacing on (a) temperature gradient and (b) CR at 650 mm/s and 0.12 mm.

Figure 33 .
Figure 33.Temperature history at constant energy density 3.67 J/mm 2 with different scan speeds and hatch spacings.

Figure 34 .
Figure 34.Surface RS at a constant energy density of 3.67 J/mm 2 .

Figure 35 .
Figure 35.Effect of the same energy density of 3.67 J/mm 2 on (a) temperature gradient and (b) CR.

Figure 36 .
Figure 36.Temperature history at a constant energy density of 2.82 J/mm 2 with different laser powers and scan speeds.

Figure 37 .
Figure 37. Surface RS at a constant energy density of 2.82 J/mm 2 .

Figure 38 .
Figure 38.Effect of the same energy density of 2.82 J/mm 2 on (a) temperature gradient and (b) CR.

4. 3 .
Part Scale Model 4.3.1.Tensile Test Tensile testing of as-built condition showcased the sample undergoing ductile failure as presented by the true stress-strain curve shown in Figure 40.The yield strength at 0.2%elongation was 650 MPa, with a Young's modulus of 171 GPa, and the ultimate strength was 1191 MPa.By fitting the data to the Ramberg-Osgood equation, the strength coefficient was found to be 1618 MPa, and the strain hardening exponent was 0.243.

Figure 42 .
Figure 42.Evolution of the part scale Mises stress at (a) 667 s, (b) 1714 s, and (c) 2310 s (d) after separation from the build plate at 220 W, 650 mm/s, and 0.12 mm.

Figure 45 .
Figure 45.Effect of hatch spacing on in-depth RS at 650 mm/s and (a) 220 W and (b) 140 W.

Figure 46 .
Figure 46.Effect of laser power on (a) in-depth RS and (b) deflection in the Z direction.

Figure 47 .
Figure 47.Effect of scan speed on (a) in-depth RS and (b) deflection in the Z direction.

Table 3 .
Process parameters for LPBF of IN625.