An efficient numerical model for predicting residual stress and strain in parts manufactured by laser powder bed fusion

Computational modeling of additively manufactured structures plays an increasingly important role in product design and optimization. For laser powder bed fusion processes, the accurate modeling of stress and distortion requires large amount of computational cost due to very localized heat input and evolving complex geometries. The current study takes advantage of a graphics processing unit accelerated explicit finite element analysis code and approximated heat conduction analysis to predict the macroscopic thermo-mechanical behavior in laser selective melting. Adjacent layers and tracks were lumped to reduce the number of time steps and elements in the finite element model. The effects of track and layer grouping on prediction accuracy and solution efficiency are investigated to provide a guidance for a cost-effective simulation. Thin-wall builds from Inconel alloy 625 (IN625) powders were simulated by applying the developed modeling approach to get the detailed residual stress and distortion at a computational speed 50 times higher than conventional approach. Under repeated heating and cooling cycles, a high tensile stress was produced near surfaces of a build due to a larger shrinkage on surface than that in central area. It is also shown that horizontal stresses concentrate near the root and top layers of the IN625 build. The predicted residual elastic strain distribution was validated by the experimental measurement using x-ray synchrotron diffraction.


Introduction
Additive manufacturing (AM) is an innovative fabrication technology to build up a part layer by layer which offers a great flexibility in product design. Among various AM methods, laser powder bed fusion (LPBF) has the advantages of a high-quality surface finish, high process reliability and ability to create complex structures [1]. Extensive studies have been focused on the porosity formation, fluid flow in the molten pool, and characterization of the powder motion, and interactions [2][3][4][5][6][7]. Experimental works [8][9][10] have found that in situ part distortion and thermal stresses can result in build failures in the form of interference, crack initiation, and so on. From the durability perspective, the AM parts should have minimized residual stresses, which is critical to maximize fatigue lives. AM by LPBF involves high costs due to expensive powder materials and relatively long building time, thus multiple trial runs should be avoided, especially for large-scale components. To date, in situ observation and management of part distortion and residual stress are difficult due to powder spreading, melting, and limited resolution of measurement techniques. Therefore, numerical models with experiment validations can be utilized to predict those quantities and provide valuable insights into process optimization.
While the numerical algorithms for metal fusion, fluid flow and solidification are advancing significantly [11], the ab initio simulation of component-level LPBF is very costly even using modern high-performance computers. The powder used in the LPBF is in order of sub-millimeter, as a result, fabricating a 25.4 mm (1 inch) cube will need 100 million powder particles and over 1000 m of laser scanning if cross sectional area of molten pool is 0.01 mm 2 . In addition, the heat source is confined to a very small region (tens of micrometer) accompanied with extremely high heating and cooling rates. The laser scanning speed is usually very high (∼1 m s −1 ) compared with direct energy deposition such wire and arc AM. The time increment to capture the steep temperature gradient is about 10 µs. Due to such multiple time and length scales, direct modeling of each powder in a part is extremely challenging in terms of computational cost and time. The reported studies focus mostly on the heat flow part rather than the stress analysis except those studies on very small dimensions such as single-track build. King et al [5] summarized the material and computational challenges in the LPBF, pointing out there are a large number of processing parameters and a broad range of time and length scales. The power absorptivity and powder interaction are critical to the heat conduction and layer formation. Ahmed Hussein et al [12] utilized the ANSYS software package to simulate the transient temperature and stress during single-layer laser scanning. They found the length of melt pool increases with the scanning speed, but the width and depth of the melt pool decreases. Matsumoto et al [13] simulated the single layer selective melting using a 2D plane stress model with consideration of temperature dependence of material properties. They assumed the powder bed and the molten area to be continuous and their Young's modulus values are very small. Wang et al [14] developed a coupled thermo-mechanical model to simulate laser cladding induced residual stresses in functionally graded material (FGM) layers. Dual-heat source with Gaussian distribution in longitudinal and transverse directions was employed based on the morphology of FGM layers. In another study, Shi et al [15] found that FGM layers can be used to reduce residual stress in a brazed joint. The location of cracking has good correlation with numerical modeling.
To predict the thermal distortion in LPBF at component level, numerical models with appropriate simplification have been proposed. The inherent strain or eigen strain method initially derived from computational welding mechanics [16] has also been used for accelerating computation at part scale. The detailed strain distribution in a local model is firstly calculated using a thermo-mechanical approach and then mapped to the global model. Keller and Ploshikhin [17] proposed a multi-scale model and employed the concept of inherent strain for simulating part distortion of powder based AM. Computation time has been greatly reduced from days to hours. Megahed et al [18] reviewed the numerical approaches at different length scales on the thermo-mechanical modeling of AM, including those with consideration of multi-scale effects and coupled multi-physics. Dunbar et al [19] analyzed the mechanical behavior of a cylindrical wall by LPBF using a code by Project Pan. The displacement profiles along the wall height were measured and a reasonable agreement between simulations and experimental measurements was obtained. In most of the relevant studies, residual stresses were not discussed because the inherent strain computed in a small model may not be suitable for stress prediction in real parts where restraint can be different. Similarly, lumped layers or peak temperature approach can be applied to save considerable computational cost. Prabhakar et al [20] applied stationary heat flux to specified layers over a time period in the numerical model. Yang et al [21] compared the line heating and layer heating model in modeling of thermal distortion of LPBF parts. The layer heating method based on Abaqus is found to be promising in thermo-mechanical simulation. Distortion results were validated by experiment whereas residual stresses were only discussed from simulation results perspective. Cao et al [22] investigated a coupled thermo-mechanical model based on Ansys for laser melting deposition (LMD) process. They found that calculation time can be saved by 85% when six layers are lumped to an equivalent numerical layer. Commercial software such MSC-Marc and Autodesk-Netfab also has the function of combining the thermal load from a heat source scale model and applying to a mesoscale hatch layer for thermo-mechanical analysis. Abaqus-AM uses the partial integration and homogenization modeling techniques to allow increased element size. Xie et al [23] and Ding et al [24] took advantage of the steady-state temperature distribution to reduce the computational time. The part distortions were generally predicted with a good accuracy, but the residual stresses were overestimated in wire and arc AM.
In the framework of transient thermo-mechanical finite element analysis (FEA), computational cost can be reduced by adaptive mesh refinement [25,26] and iterative substructure method [27,28] without compromising solution accuracy. More recently, Ma and Narasaki [29] proposed a hybrid explicit/implicit FEM for welding stress simulation and hybrid solid/shell modeling to save computation cost. For the first step of stress analysis, adaptive mesh and dual-mesh methods [30] can be employed to get the temperature history. To accelerate the heat conduction analysis, Luo and Zhao [31] developed a mesh refining and coarsening scheme to solve finite element equations for temperature field. With optimized multiprocessor parallelizing, the layer and part level powder selective melting can be simulated 12-18 times faster compared with traditional simulation methods. Gan et al [32] conducted a comparative study on the simulation models for temperature field and cooling rates. They concluded that the melt pool size predicted by heat conduction model, thermal-fluid model, and thermal-fluid vaporization model all have shown reasonable agreements with experimental results [33], with only difference for cooling rate which determines the solidification process. Lee and Zhang [34] proposed a numerical approaching by combining a powder packing model with a transient heat and fluid flow model. The resulting temperature and solidification rates were employed to predict the grain size and molten pool. As mentioned earlier, the solution domain is usually quite limited due to the tremendous computational cost when multi-physics are involved. Promoppatum and Yao [35] studied the effect of scanning length and energy input on residual stress in Ti-6Al-4V AM. The computational domain in 2.5 mm cube with 20 scanning tracks are explored.
In the present study, we investigated the residual stress and strain in relatively larger parts fabricated by LPBF. The dynamic laser scanning passes are organized into groups and progressively activated to simulate the residual stress and distortion of LPBF-fabricated parts. The effect of track and layer grouping scheme on the prediction efficiency and accuracy was quantitatively evaluated by a numerical model. To simplify the heat flow modeling, an analytical solution of temperature field under a fast-moving heat source was employed as temperature input for stress analysis. The heat efficiency of the grouped laser scanning was calibrated to ensure an equivalent volume of molten pool. Recent developments have been made to use explicit FEM as the stress analysis solver and accelerate the computation based on graphics processing unit (GPU) techniques [36][37][38]. An in-house explicit FEM code [38] based on GPU acceleration was employed to model the residual stress and strain distribution in two LPBF builds. Residual elastic strain modeled by developed simulation approach was compared with experimental measurement by x-ray diffraction.

Thin wall Inconel alloy 625 (IN625) build for residual stress study
The Air Force Research Laboratory (AFRL) released a competition program 'AM modeling challenge series: macroscale process-to-structure' in August 2019 to test and validate the modeling capabilities of LPBF from academy, businesses and institutes both in United States and internationally [39][40][41]. Multiple IN625 builds were deposited with the LPBF process using an EOS M280 metal 3D printer. Gas atomized IN625 powder was used as the LPBF stock. The particle size of this IN625 powder is between 6 and 60 µm with a bimodal distribution. The plain carbon steel substrate (250 mm × 250 mm × 30 mm) for deposition was preheated to 80 • C before printing. The laser spot diameter is 0.1 mm with a Gaussian distribution. The used LPBF parameters include a laser power of 300 W, a scanning speed of 1230 mm s −1 , and a hatch spacing of 100 µm. The laser scanning follows the rastering path, and the stripe length is 10 mm. The scanning vector is rotated 67 • after completion of each layer. The layer thickness was 40 µm. No post build heat treatment was applied to the builds after deposition. The calibration wall builds attached with the steel base were machined by electrical discharge machined for the residual stress measurements. In this study, two builds A46 and A51 were selected for FEA simulation. The dimensions of the two builds are illustrated in figure 1(a), A46 (10 mm × 1 mm × 25 mm) and A51 (50 mm × 1 mm × 25 mm). The actual A46 wall sample for calibration (averaged thickness 0.8 mm) is shown in figure 1(b).

Measured material properties and elastic strain map
Mechanical properties of the IN625 builds are essential parameters of residual stress and distortion simulations. The yield strength and tensile strength of the IN625 builds in different build orientations were tested by General Electric (GE) Global Research Center sponsored by AFRL [42]. Figure 2(a) shows room temperature tensile properties of the as-built IN625 part tested in different orientations (X, Y, Z). The build direction (Z-direction) shows the lowest strength at room temperature comparing with the other two orientations. The tensile properties of the as-built parts in the build direction (Z-direction) were also measured at different temperatures (25 • C, 538 • C, 816 • C, 982 • C, and 1093 • C), as shown in figure 2(a). The measured elastic modulus and tensile strength of the AM IN625 builds in Z direction as a function of temperature is shown in figure 2 In the validation dataset, residual elastic strains within the two selected IN625 builds were provided based on experimental measurements. The in-plane strain components horizontal elastic strain (e xx ) and vertical elastic strain (e zz ) in the builds were measured by the energy dispersive diffraction (EDD) technique using x-ray synchrotron at the advanced photon source at Argonne National Laboratory [40].

Finite element model
The one-way coupled thermo-mechanical FEA is generally used to analyze stress and distortion in AM. The computational procedure can be found in literatures [24,43,44]. In the present study, heat conduction analysis is approximated by an analytical solution of a point heat source considering the fast scanning speed in LPBF process. The molten pool shape and residual stress prediction accuracy will be employed to justify the suitability of this simplification for stress analysis. The elastic modulus and tensile strength of AMed IN625 parts are different in Z-direction and X and Y directions (figure 2). To achieve higher accuracy prediction in stress and strain, the orthotropic plasticity hardening model [43] was employed in the modeling of material behavior of LPBF build. For simplicity, the same yield stresses are assumed in the horizontal x-and y-directions, and the stress ratio between vertical z-direction yield stress and horizontal yield stress is 0.868 based on measurement data at room temperature. Other material properties such as thermal expansion coefficient and Poisson's ratio can be found in [43].

Heat input calibration
The temperature history resulting from the laser heat input is important to the stress evolution in metal AM. One of the key factors is net heat input in the powders considering the multiple reflection of laser beam and complex interaction among powders. According to the study by Boley et al [45], Zanzarin [46], the laser absorptivity is 0.67 ∼ 0.75 for iron particles if the sizes of powder have a uniform distribution. Here, the net heat input is calibrated as 0.42 considering the measured density ratio (0.6) of compact powder with respect to deposited metal. The average porosity rate calibrated from 2D cross section of sample A46 and A51 are 0.18% and 0.29%, respectively. A sensitivity study has shown that the small porosity rate (<1%) will have very limited effect on residual stress prediction accuracy. Thus, the powders in the vicinity of laser radiation are assumed to become continuum solid after melting in the numerical model. For simplicity, powder sintering process is not simulated, and the analytical solution of point heat source on the surface of semi-infinite body is employed to approximate the temperature distribution of laser selective melting as shown in equation (1). For fast moving heat source, it would be more convenient to use an approximated solution of the Rosenthal equation by assuming heat spreads in transverse direction only [47] as described by equation (2). where q is the net laser power, λ is the heat conductivity, R is the distance from temperature evaluating point to heat source, x is the vector component of R along heat source moving direction, r is the vector component of R in the direction perpendicular to the laser scanning track, a is the thermal diffusivity, v is the moving speed of heat source and t is the traveling time of the heat source relative to an evaluating point. The transient temperature at a point can be easily superposed for multiple scanning around the point, which is the case for LPBF.
Tang et al [48] compared the predicted melt-pool widths of single beads by the Rosenthal equation with experimentally measured data by Yadroitsev et al [49]. They found a good agreement in width but 30% difference in cross sectional area, especially for lower power density cases. In this study, the heat conduction in single-track laser scanning was simulated by FEA by JWRIAN developed in Osaka University. The temperature-dependent material properties [43] were used for heat conduction analysis considering the thermal resistance of powders [50]. The Marangoni effect was indirectly considered in the thermal model to calibrate heat input. Due to strong convection induced by Marangoni flow, an enhanced heat conductivity by three times was applied for the values beyond the melting point in the FEA modeling of heat conduction. The deposit metal is 0.15 mm in width and 0.5 mm in height. The adjacent powder area has equivalent thermal conductivity at condition of minimum contact of powder 0.28%. For analytical solution, the thermal physical properties such as heat conductivity and specific heat was chosen at 1000 • C as baseline. By fixing thermal diffusivity, the thermal conductivity was adjusted to examine the melt pool variation. As shown in figure 4, the heat conduction as predicted by the Rosenthal equation is compared with FEA. At equivalent thermal conductivity 0.022 W m −1• C −1 , the analytical melt pool length is close to that by FEA prediction. On the other hand, a significant amount of the latent heat of fusion will be absorbed/released during melting/crystallization. The single-track model was used to test the effect of latent heat on temperature distribution and residual stress evolution. It was found that the melt pool in the case with latent heat is about 9% longer than that in model without latent heat, because the latent heat is released during rapid solidification of melt pool. Nevertheless, the stress results confirm that temperature difference induced by latent heat does not have obvious impact on residual stress. For this reason, the latent heat is neglected in the simplified heat conduction model for stress analysis.
To validate with available experimental data, a NIST benchmark example [32,33] of single bead-on-plate laser weld on IN625 is used for reference as shown in table 1. The heat input has the same energy density 0.24 J mm −1 (energy per length) under the condition of laser power 195 W and travel speed 800 mm s −1 . The analytical solution gives a reasonable prediction on the molten pool size for the single bead laser scanning. It should be noted that the thermal conductivity has large dependence on temperature which indicates the temperature at region away from the heat source may be underestimated. The difference in temperature is partly compensated by the assumption that no radiation and convection is present in the thermal model. Despite such simplification, the molten pool size has a good agreement between the analytical and numerical cases.

Grouping of track and layer (GTL)
To predict distortion and residual stress at a component level, both length scale and time scale can be accelerated to a reasonable extent without comprising accuracy substantially. One of the options would be lumping multiple scanning tracks and layers, in which case, mesh size and time increment can be greatly increased. While this computational model may be less efficient than inherent strain or eigen-strain based models, more reliable prediction in the residual stress is anticipated since the elastoplastic response of material is solved. To use the Rosenthal solution in a laser scanning with a back-forth pattern, the approximated solution for a high-speed welding was employed. As shown in equation (2), the transient   temperature at a location is related to the laser scanning path location and heat source traveling time. The temperature at a material point is determined from the peak temperature induced by each individual heat sources rather than the summation of temperatures induced by those heat sources. In other words, there is no thermal interaction between each track or layer that are lumped together for stress analysis. Moreover, the effect of cooling time on temperature gradient is considered for each layer. An increased cooling time is added for the lower scanning paths based on the layer grouping factor G L . The effect of grouping scheme on solution accuracy and efficiency is investigated to find out the optimum grouping parameters. Here, a sensitivity study is performed on six simulation cases of grouping scheme as shown in figure 5. A 2.0 mm × 0.5 mm × 0.48 mm box model is simulated with laser scanning  Figure 6 shows the temperature field by proposed analytical solution at the last three layers in case A. The molten pool shape predicted by the analytical solution is similar to those by the FEA result and the original Rosenthal solution. The melt pool extends in the reverse direction as the laser scanning continue to the adjacent track. The temperature field is smooth and continuous in each layer scanning with different orientations, which is important to ensure the convergence and accuracy of stress analysis. After scanning of each layer, a cooling time is set to allow the peak temperature of print parts cool down to below 80 • C as specified in the physical experiments. Figure 7 shows the residual stresses calculated by the proposed method with different grouping schemes. The horizontal stress (stress component along X direction) exhibits high tensile stress near the top surface and compressive stress near the bottom layer, which is typical in longitudinal welding of beam structures. The vertical stress (stress along Z direction) has high stress concentration at the lower part due to the strong constraint applied at the bottom. From case A to case E, the resolution of stress becomes lower as the mesh size and grouping factor increase. Nevertheless, the key features such as stress pattern and magnitude are maintained to a satisfactory extent, all showing high tensile stress in horizontal direction near the top layers and low compressive stress in vertical direction at the bottom. To quantify the effect of mesh coarsening and grouping scheme on stress accuracy, the simulation results along the center line on the top layer of the printed block are compared in figure 8.

Computational accuracy and efficiency by GTL approach
From the stress profile in figure 8(a), the baseline model without grouping demonstrates a relatively constant stress magnitude (800 MPa) except the location near two ends. A minor wavy pattern can be observed due to the track-by-track laser scanning process, reducing residual stress in adjacent tracks by annealing effect. Such phenomenon was also observed in other studies [25] using the similar powder material (Inconel 718). With track and layer grouped, the stress magnitude has larger oscillation because the adjacent two or three layers are heated at the same time using the same scanning pattern, and deviation from baseline case becomes larger especially when grouping factor exceeds three. The horizontal displacement profile in figure 8(b) has similar phenomena of wavy pattern as the grouping factor increases, but the shrinkage measured from the two ends coincides well among the four cases. The overall horizontal shrinkage of the build predicted by the GTL model has less than 10% difference from the baseline case (case A).
The maximum stress and in-plane shrinkage obtained from different cases are summarized in table 2. The grouping scheme and factor have an obvious influence on the computational time and solution accuracy. With only scanning tracks lumped (cases B and C), the value and location of peak stresses are both well reproduced compared with the baseline model. In those cases, the computation time is reduced to nearly a fraction (1/G T ) of the cost by case A. With the laser scanning lumped in two directions (cases D and E), the mesh size can be increased because the dimensions of molten pool and heat affected zone are larger. Increasing the mesh size from 0.02 mm (case A) to 0.06 mm (case E), the scale of finite element model is drastically decreased from 50 000 to 2500 which amounts to a reduction factor of 3 3 = 27. Accordingly, the computation time can be reduced by a similar factor if a fixed time increment is used which is the case for most implicit FEMs. For the explicit FEM, this will bring in an additional factor of 3× since the stable time increment is also increased by mesh size. In the simulation of a real-scale model, grouping scheme illustrated by case E was employed to balance the computational time and accuracy. Figure 9 shows the predicted elastic strain and residual stress after the fabrication of A46 build. The contours of stress and strain are basically symmetrical about the central line due to the rectangular shape of build and back-and-forth laser scanning pattern. Comparing the stress contour with the strain contour, a similar distribution can be observed, indicating the vertical strain has a dominating effect on the vertical stress component. The maximum strain near vertical edges of thin wall is about 0.5% which exceeds the elastic limit of IN625, thus stress is expected be high (>700 MPa) as a result of strain hardening. The residual elastic strain and stress in vertical direction has significant variations in horizontal plane. The vertical stress is highly compressive in the center of the LPBF build but highly tensile at corners. Such stress pattern is very consistent with those in L-shaped rectangular specimen measured by combined digital image correlation method and neutron diffraction method [5]. The stress concentration along the vertical build direction is harmful to the material properties and the structure integrity of prints, which may partly explain the lowered yield strength and properties in Z direction. The internal residual strain appears to be compressive as a result of stress balance in the vertical direction. Figure 10 shows the residual elastic strain predicted by the numerical model and those measured by x-ray synchrotron. Note that, the elastic strain is the averaged quantity through the thickness according to nature of the characterization approach and interaction volume described in section 2.2. There is high tensile elastic strain (horizontal component) on the top layer due to the fast heating and cooling of the laser scanning. Intermediate part of the build has a low elastic strain because repeated laser scanning thermal cycles have certain extent of heat treatment effect on the earlier printed layers. The bottom layer retains the high residual stress and strain since the restraint from the substrate is large which prevents the shrinkage to occur. On the other hand, the vertical elastic strain exhibits a high tensile strain near the vertical edges covering almost all layers of the build. The strain distribution by x-ray measurement is not as smooth as that by prediction near the central area of the build probably due to the rotating pattern of laser scanning and simplification of the model. Nevertheless, the model predicts the pattern and magnitude of elastic strain reasonably well. For more quantitative comparison, the strain profile along three lines are plotted in figure 11. The horizontal strain is smaller than vertical strain especially on the edges (lines 1 and 3), resulting in higher residual stresses in vertical direction. From results along line 2, the two strain components are both small (∼0.1%) through most printed layers except the top and bottom part. Compared with experimental data, the model seems to underestimate the horizontal strain in the center and overestimate that on the edge. Such difference in results comparison can be attributed to the simplified laser scanning pattern and material anisotropic properties (Poisson's ratio and Young's modulus). In general, the horizontal strain is small compared to the vertical strain component and elastic limit.

Dependence of residual stresses on part dimension
The two builds A46 and A51 with different dimensions were simulated by the proposed computational model and compared with the available x-ray synchrotron measurement data. The height of the two models are  both 25 mm, while A51 is wider and thicker (30 mm × 5 mm) than A46 (10 mm × 1 mm). The same laser scanning schemes and parameters were used for the builds. For specimen A51, the vertical strain predicted by the developed model is plotted in figure 12, along with strain contour on two sections in the center of build. Similar to the strain map of A46 (figure 9), the internal region exhibits highly compressive strain which is close to elastic limit. With the increased width and thickness, all regions near surface have high tensile strain other than the corner region in the case of A46 build. The restraint in both x and y directions are high in the case of A51 build, whereas the A46 build is closer to plane stress state due to its much smaller thickness (1 mm). The strain gradient through thickness in large, caused by the restraint variance between surface and interior of the build. The vertical strain exhibits a similar tensile-compressive-tensile profile across the lateral x and y directions, with maximum strain about 0.5% on the edges. The residual stress on the top layer and bottom layer are plotted in figure 13 for the two builds with different dimensions. The horizontal stress on the top layer has almost the same type of relatively uniform distribution and peak stress magnitude (750 MPa) except at the two ends. This can be explained from the welding mechanics perspective by regarding the LPBF of top layer as welding on the edge of a specimen. After experiencing rapid heating and cooling, the material of top layer tends to shrink under restraint by existing build beneath the layer. The same level of heat input will result in similar level of inherent strain as well as horizontal stresses. The stress pattern on the bottom layer shows higher stress near the surface but lower stress at the center. The difference from that in top layer is mainly caused by the high-level vertical elastic strain near the surfaces. The stress distribution is 2D on the top layer whereas a 3D stress space is expected at the bottom layer. The peak horizontal stresses in A51 are higher than those in A46 despite of certain similarity in stress distribution. This is probably due to the enhanced restraint in y and z directions by thicker specimen that increases the size of region with highly tensile residual stress in vertical direction. From the subfigure at bottom, the vertical stress with a high magnitude has been produced in the corner regions of the two specimens. In the case of A51, such high stress region extends to nearly all surfaces, changing the stress distribution in horizontal direction.
The comparison of residual elastic strain between model and experiment is shown in figure 14. Considering the nearly symmetric distribution of elastic strains, results on half part of model are plotted to facilitate direct comparison in the same figure. The horizontal strains predicted by model and measured by x-ray both feature compressive strain on two sides of the build and tensile strain on top and bottom layers. Peak horizontal stresses appear on the top and bottom layers as demonstrated in figure 13. Inside the build, the averaged elastic strain ranges from about −0.1% to 0.1%, indicating the stress is relatively lower than that in surface region. The vertical elastic strain is in highly tensile state on two sides of the build, and the hot spot region is about 2 mm wide at the root of the build. The averaged strain is compressive in the internal region of the build, which comprises a very thin layer of tensile strain as confirmed by the model previously. It is also interesting to note that, the map of horizontal strain inversely correlated to that of vertical strain, where the region with high horizontal strain has low vertical strain and versa vice. Overall, the model well reassembles the residual elastic strain distribution in experimental specimens, more accurate prediction would require consideration of actual laser scanning pattern, fluid-flow and mass transfer, powder dynamics and improved heat conduction modeling.
The analyses by developed in-house explicit FEM code are completed within a short time frame. For A46 and A51 simulation models, the computational time is 2 and 26 h, respectively. All computations were performed on a workstation equipped with Intel Xeon Silver 4110 Processor (8 Core) and NVIDIA Tesla Pascal 100 GPU. Due to large number of elements and time steps involved LPBF model, such simulation would require weeks to months of computation effort by commercial software. Ren et al [51] simulated a deposition in a domain 72 mm × 7 mm × 0.5 mm (single layer) which has comparable build volume of A46 (25 mm × 10 mm × 1.0 mm). With similar mesh size used in this study, the total simulation time of thermomechanical analysis used by COMSOL software was more than 105 h. In contrast, the simulation time of the modeling case A46 was only 2 h which is less than 1/50 of Ren's numerical example. Adaptive mesh method and multi-GPU technique can be applied to further accelerate the computation, providing a cost-effective tool to evaluate the residual stresses in LPBF manufactured parts.

Concluding remarks
A multiple track and layer grouping approach was developed to accelerate the FEA simulation of laser selective melting. The temperature field was approximated by an analytical solution using the Rosenthal's equation, which predicts the molten pool with a reasonable accuracy. The computational time is greatly reduced by the GTL technique and GPU acceleration. The following conclusions can be drawn based on this study.
(a) For laser selective melting, multiple tracks and layers can be grouped together with temperature history determined by each fast-moving heat source. The optimum number of tracks and layers are two-three from the computational accuracy and efficiency perspectives. (b) The length and width of molten pool predicted by the simplified analytical solution are both about 18% lower than prediction by FEA and experiments. In addition, the calculated depth of molten pool is 15% smaller than FEA prediction but 30% larger than experimental data, indicating a reasonable prediction of the temperature field for fast moving heat source.