Introduction

Additive manufacturing (AM) has been of interest from as early as the 1980s [1]. The process requires that the desired shape be divided into relatively thin layers (or slices) and added consecutively on top of each other. Initially used to produce rapid prototypes, significant research and development over the past decade means that AM can now be used as a near net shape method of manufacture. The technique boasts a wide range of applications from weight-reduced aerospace parts [2] to bespoke artificial joints [3]. Between 1990 and 2015, the growth rate of the worldwide revenues generated from all AM products or services was 25.4% [4]. This increase in interest was due to the wide range of benefits the technique has to offer. AM provides methods of waste reduction [5], customization opportunities [6], and part repair [7]. Furthermore, it eradicates the need for the creation of geometry specific dies and provides new methods of creating even more complex geometries.

One particular AM process that offers the capability of the production of complex structures is laser powder bed fusion (LPBF), which is also commonly referred to as selective laser melting (SLM). SLM is a powder bed fusion process that uses a laser as a heat source as opposed to an electron beam or wire arc. The use of a powder bed enables the fabrication of lattice structures; features that could not be manufactured by traditional methods such as casting. SLM is one technique that can also be used to reduce the weight-to-strength ratio of a part by manufacturing topology optimized parts. However, due to the relative infancy of the technology, there are still a number of challenges with the process that are preventing its widespread introduction into industry.

The procedure has a number of limitations, some of which challenge the reliability of the build parts and others that constrain the manufacture process itself. For example, certain AM procedures such as SLM are limited by the size of the build chamber. These processes can also exhibit very long manufacture times, due to the slow scan speeds and small layer heights. Other issues that can be more widely applied to all AM procedures include surface roughness, lack of fusion, and distortion.

Primarily, the main challenge with all AM procedures is the reliability and repeatability of parts [8, 9]. Due to the thermal cycles that are experienced by the part during the build process, residual stresses are induced, resulting in part distortions. In addition, the thermal field induces changes within the microstructure of the material, resulting in altered material properties. These factors cause unreliability in additively manufactured parts. Alterations in mechanical properties and unknown residual stresses mean that the suitability of the part is unknown until tested. Therefore, non-destructive testing methods are required to determine the material properties of the built part and any possible defects, such as pores caused by lack of fusion; however, the lack of inspectability of complex parts by means other than computed tomography limits scalability.

The key to unlocking the potential of this new technology is understanding every aspect of the process and its effect on the material behavior, so that parts can be manufactured in a fit-for-purpose state. Numerical modeling plays a large part in enabling this by providing a method of testing the effect of physical parameters without the cost of experimental trials. Once a validated model that can predict the impact of these parameters has been established, it can be used to investigate of the tailoring process parameters to achieve desirable material properties and minimize distortion and residual stress.

Previous work on additive manufacturing has included the research of modeling methods to predict distortion and residual stresses [10,11,12,13,14], melt pool geometry and energy penetration depth [15, 16], microstructure [17, 18], or to determine the effect of different process parameters [19, 20]. Due to the complexity of the process, simulation techniques vary a great deal, depending on which particular AM process is being implemented and what feature or parameter is the central interest of the study. A relatively detailed review of existing modeling approaches was recently summarized by Luo and Zhao [21].

In this study, the focus is on a method of prediction of residual stresses for the selective laser melting process. In order to model residual stresses, first a sufficiently representative thermal model is required.

A number of heat source models have been used to model the laser. Song et al. [22] undertook an investigation into residual stresses and microstructure using powder bed methods with a nickel superalloy. A Goldak double ellipsoidal model was used to represent the moving heat source. Meanwhile, Hussein et al. [23] use the more common Gaussian model in their analysis of temperature and stress fields. While Hussein et al. [23] use only a single layer in their model, simulations which include multiple layers must exhibit a method of material deposition. Song et al. [22] use the element birth and death technique, in which sets of elements are activated/deactivated within new steps. This method is implemented by many others including Fu et al. [24]. Alternative methods include the quiet element method demonstrated by Michaleris [25].

Material properties, specified within the analysis, should also be temperature and state dependent [26]. For example, Huang et al. [27] establish an effective thermal conductivity for powder particles, and also use a linear mixing rule to take porosity of the powder into account within its density.

When considering the simulation of residual strains, a number of approaches are available. Li et al. [11, 12] have developed a multi-scale model. Multi-scale models have become of increasing interest due to their computational efficiency. The model shows the implementation of a Gaussian heat distribution on a micro, meso, and macro scale. The benefit of this method is that it allows for the implementation of specific scan strategies. Other approaches include an analytical method, demonstrated by Fergani et al. [10]. Furthermore, Wu et al. [13] implement a sequentially coupled thermo-mechanical analysis through Abaqus, similar to the work presented here. Sequentially coupled analyses save a significant amount of computational time [26] without sacrificing a reasonably accurate representation of the underlying physics. Denlinger et al. [28] also use a sequentially coupled analysis in their model of electron beam direct manufacture and again in reference [29]. In this analysis, the weak coupling was used as the plastic strain energy was significantly small in comparison to laser energy. More recently, Williams et al. [14] have also used a sequentially coupled analysis through Abaqus, using a ‘block dump’ approach for material deposition that will be used within this work. This is an efficient tool that also reduces computational cost of the analysis.

In this paper, a specific approach for simulating residual strains within an Inconel 625 part created by SLM is described. This was undertaken as part of the AM Benchmark 2018, challenge 1 to (AMB201801). In section “The Challenge”, an outline of the challenge is provided, followed by a description of the modeling approach in section “Modeling Approach”. Finally, a discussion of the results is presented in section “Results and Discussion”.

The Challenge

For this section, note that all details regarding the design and setup of the experiment have been sourced from reference [30], along with some independent correspondence with the challenge coordinators.

The part geometry for this challenge was a simple single cantilever structure of dimension 75 mm, 5 mm, and 12.5 mm in length, width, and height, respectively. A depiction of the geometry can be seen in Fig. 1.

Fig. 1
figure 1

The full geometry of the build part as given by NIST [30]

For each build process, four parts were built on the same substrate. A full schematic of the layout and dimensions can be found in reference [30]. However, as spacing is sufficient between parts such that the build of an adjacent part has a negligible effect on the thermal history of the current part, the model is reduced to a single part, saving on computation time. Nevertheless, it is important to consider the entire build layout when calculating the cooling time between layers, as this will have a first-order impact on the thermal history of the part and consequently the residual strain field. The reduced model is a single part, mounted on a section of the substrate. The section is taken to be 81 mm × 12.7 mm × 11 mm (length × height × width).

Within the benchmark, the part was built in IN625 and 15-5PH Stainless Steel as part of the experimental tests held at the NIST laboratories. The substrate in both cases was the same alloy as the build part. In this study, IN625 was chosen as the build material. The material properties used within the model can be found in section “Material Properties”.

The manufacture of parts, within the benchmark, was undertaken on two different machines: an NIST-built machine and an EOS M270. It is assumed from here on out in this analysis that the EOS M270 is used. The exact build parameters used during the process can be seen in Table 1.

Table 1 Build parameters [30]

The scan strategy used is outlined in great detail within the challenge description, including laser on and off times. Firstly, each layer begins with a contour scan, creating only the outline of a single part. An infill scan of the current part then occurs, before moving on to the next part within the build. From Table 1, it can be seen that the infill and contour scans have different laser properties. In addition, the infill scan strategy alternates, between horizontal and vertical scans, depending on whether an odd or even layer is being completed. This is a common technique used within AM in order to reduce lack of fusion within the build.

The parts described above were built to these specific parameters at the NIST laboratories. This was in preparation for experimental testing as comparative data for simulations.

Multiple tests were undertaken on the completed build part and many features observed, including distortion, residual strain, and microstructure properties. This paper focuses solely on the prediction of the residual strain measurements.

The evaluation of residual strain was undertaken using X-ray diffraction (XRD) measurements across the central plane of the build. Additionally, measurements were only taken along the thickest legs. The measurement plane can be seen in Fig. 2.

Fig. 2
figure 2

Plane for XRD measurement [30]

For additional information on the experimental measurement procedure, please refer to reference [31].

Modeling Approach

In order to achieve the prediction of residual strains, a numerical simulation was undertaken using the additive manufacturing process Apps from DASSAULT SYSTEMES [22]. Details of the modeling approach in the software package can be found below in section “Software”. The software offers several modeling approaches including thermo-mechanical and eigenstrain/inherent strain methods. In this particular work, the authors have adopted the sequentially coupled thermo-mechanical procedure.

Sequentially coupled analyses rely on the fact that the first analysis is independent of the second. Hence, in this particular case, the thermal history of the build can be assumed to be independent of its mechanical response [29]. This form of analysis is preferential to a fully coupled analysis as it saves a considerable amount of computational time and cost. Sections “Thermal Simulation” and “Mechanical Simulation” outline the approach to each part of the analysis.

Software

Initially, the DELMIA Powder Bed Fabrication App was used to generate the slicing, recoating, and laser trajectories of the part. This was done to comply with the specifications of the challenge, including laser speed, power, and infill pattern. The actual analysis was completed in the SIMULIA Additive Manufacturing Scenario App, based on the Abaqus 2018 finite element solver. The app offers the ability to use either a thermo-mechanical method or an eigenstrain based method. Within this study, we used the thermo-mechanical approach.

As implemented in Abaqus 2017 [32], key features of the physics-based framework that facilitate a more representative analysis of additive manufacturing processes include (1) the ability to use finite element meshes that do not conform to the actual powder layer thickness; (2) explicit use of the spatial and temporal data associated with laser scan paths and powder recoating for defining material activation in the model and integration of the moving heat source; (3) an “intersection module” that slices the laser scan path within the finite element mesh to distribute the heat source over a given increment of time; (4) progressive element activation to track the specific volume of an element that is either completely or partially filled with material (integration points) that are active within a given time increment; and (5) progressive heat loss through convection and radiation by determining the evolving lateral and vertical free surfaces as a consequence of the layer-by-layer activation process. These specific features and their application within the current finite element model are described in more detail in subsequent sections.

Geometry and Mesh

The single cantilever used in this study can be seen earlier in Fig. 1. In order to replicate this geometry exactly within the model, the STL file was taken directly from the challenge description [30] and converted into a solid part within the Abaqus model. A second part was created in the model to represent the build plate. This is an important aspect that must be accounted for in order to achieve a reasonable thermal history prediction, due to the large heat sink effect induced by the relatively large size of the substrate in comparison to the build. The substrate dimensions were given by NIST in the task description [30]. However, to model the entire build plate would have been computationally expensive, especially as four parts were built on each build plate. Moreover, each part was thermally isolated from the other parts; that is, the spacing is sufficient such that the build of an adjacent part has a negligible effect on the thermal history of any other part. Therefore, within this model, a reduced substrate was used and only a single part was modeled. The reduced substrate was still sufficiently large (with the actual through-wall thickness modeled) to account for the heat sink effect of the substrate. The dimensions of this section can be found above in section “The Challenge”. The two parts were connected by a mesh “gluing” technique constraint within the interaction model. This ensures that the heat transfer flows from the build material through to the build plate directly.

In order to achieve a regular, structured mesh, the build part was partitioned. Throughout the thermal analysis, linear hexahedron (DC3D8) elements were used. An identical mesh was used in the structural analysis, to allow the analysis to take place on a compatible mesh (see section “Mechanical Simulation”). The characteristic element size of the build was 0.20 mm (Fig. 3). Hence, each element accounted for approximately 10 real layers. Single element layers were unfeasible due to the large number of layers within the model. The element size was chosen as compromise between accuracy and run time. The partial integration and homogenization modeling techniques in the Abaqus 2018 solver were leveraged to mitigate the discrepancy between element size and layer thickness.

Fig. 3
figure 3

Finite element mesh of part in 2 directions

Material Properties

As described previously, IN625 was the build material in this study. Due to the abundant use of this nickel superalloy within industry, material property data is widely available for this particular alloy. However, within this application to a powder bed process, temperature-dependent data is required. The large thermal gradients that take place in the manufacturing process mean that the material experiences a wide range of temperatures and physical states.

Table 2 Temperature-independent material properties

The above (Table 2) gives the values of some temperature-independent properties of the material. The temperature-dependent data used for the thermal conductivity, specific heat capacitance, Young’s modulus, Poisson’s ratio, and coefficient of thermal expansion can all be seen in Fig. 4a–e [33]. In general, for manufacturing process simulations, one accounts for the temperature-dependent plasticity behavior of the material under consideration. For example, temperature-dependent stress-strain curves would be implemented in the finite element solver to accurately account for the softening response at elevated temperatures. However, due to the range of temperatures that are experienced by the material, as a result of the time integration method used, the plastic behavior exhibits only a small dependency on temperature. Moreover, Li et al. [11, 12] and Wu et al. [13] have shown that for similar materials and geometric configurations, the temperature-dependent behavior is of second-order importance. Wen et al. [34] have also shown reasonable success predicting residual stress with temperature-independent properties. Therefore, for the simulations in this work, temperature-independent yielding behavior was implemented, with a bilinear, orthotropic plasticity hardening model as shown in Fig. 4f [35]. The material’s yield stress is 725 MPa in the horizontal x- and y-directions and 615 MPa in the vertical z-direction. The stress ratio is 0.8483. The stress corresponding to a plastic strain of 0.35 is 990 MPa in the horizontal directions, and the stress in the vertical direction is scaled with the same stress ratio as that used for the yield stress. This modeling assumption—using a temperature-independent yield behavior—has been justified by the accuracy of the predictions presented herein.

Fig. 4
figure 4

Temperature-dependent material properties and plasticity: a thermal conductivity, b specific heat, c elastic modulus, d Poisson’s ratio, e coefficient of thermal expansion, f orthotropic plasticity

Thermal Simulation

The initial step of the analysis was to simulate the thermal history of the process. This is important as it is used to drive the mechanical responses and development of residual strains within the build. In order to do so, a number of important features of the additive manufacturing process had to be modeled effectively. Firstly, the moving laser heat source needs to be modeled in such a way that is representative of the shape, the power, and the movement of the laser.

In this model, a single thermal analysis step, with duration slightly larger than the total build time of the part, is undertaken. This is to allow some extra time, for cooling, post build. The laser is modeled by a moving concentrated point heat source. This method applies the heat flux at a singular moving point, and the heat is then dissipated through the model by the underlying rules of conduction. This approach is valid when the size of the heat source (or melt pool) is small relative to the characteristic element size. In such a case, the shape of the heat source is not as important as the total magnitude of heat flux; however, it is important that the precise temporal and spatial locations of the heat flux are taken into account within the modelling approach. The energy deposition into the system is computed by taking into account the actual path of the heat source. The toolpath-mesh intersection module provides the information pertaining to the energy deposition. Thus, for any given solution increment of time duration dt, the mesh intersection calculation identifies all elements (and their associated mass properties) through which the laser passes. The concentrated point heat source is then distributed according to the mass properties of all of the intersected elements. The movement of the heat source was implemented with the use of an event series replicating exactly the movement of the scanning strategy. For more detail on this approach, please see references [36, 37] where the thermal approach is fully described. Images of the scan strategy used can be seen in Fig. 5.

Fig. 5
figure 5

Event series defining the laser movement

Another key attribute of the selective laser melting process is the material deposition. The arrival of material needs to be accounted for as it has a large impact on the heat loss effects within the system. The simulation of material deposition has been achieved in the past through the implementation of one of two main methods: element birth techniques or quiet element methods. The former of the two methods uses a model change within each new step to introduce a new set of elements. The quiet element method assigns all elements with suppressed material properties until activated by the laser.

Within this simulation, material deposition was achieved via the element birth technique within the ABAQUS AM app. This technique allows material to be added to the model within a single step. An event series was developed detailing the movement of the recoater blade and implemented within the app. The material is progressively activated using the inbuilt progressive element activation techniques, according to the recoat times given while leveraging spatial and temporal homogenization techniques mentioned above.

Further details were taken into account using boundary and initial conditions. In order to replicate the substrate heating within the model, the substrate was assigned an initial temperature of 80 °C, along with a boundary condition on the bottom surface at 80 °C to maintain this temperature. Moreover, the build part was assigned an initial temperature of 40 °C as this was taken as the chamber temperature.

Heat loss was taken into account through convection and radiation, given, respectively, by the equations below.

$$ {\displaystyle \begin{array}{l}{q}_{conv}=h\left(T-{T}_{\infty}\right)\kern3.50em \left[1\mathrm{a}\right]\\ {}{q}_{rad}=\varepsilon \sigma \left({T}^4-{T_{\infty}}^4\right)\kern2.00em \left[1\mathrm{b}\right]\\ {}\begin{array}{cc}\mathrm{Eqn}.1\ \mathbf{a}\ \mathrm{Heat}\ \mathrm{loss}\ \mathrm{due}\ \mathrm{to}\ \mathrm{convection}.\mathbf{b}\ \mathrm{Heat}\ \mathrm{loss}\ \mathrm{due}\ \mathrm{to}\ \mathrm{radiation}& \end{array}\end{array}} $$

Here, h is the heat transfer coefficient taken to be 18 W/m2 K [29] and ε is the emissivity taken to be 0.45. T is the ambient chamber temperature and σ is the Stefan-Boltzman constant. These conditions were assigned to the evolving lateral and vertical free surfaces (calculated based on the specific activated elements for a given solution increment). Heat loss through the evolving lateral surfaces to the surrounding powder was assigned the same convective heat transfer coefficient as the evolving top surface of the part. In preparing the model, a number of sensitivity studies were undertaken and it was observed that the assumptions with regards to heat loss to the surrounding powder (which in reality would be lower than the top surface) only weakly influenced the results. This is likely a consequence of the thermal model analyzing relatively large solution increments. A thermal model focusing on melt pool size predictions would require smaller time incrementation and therefore would experience a larger influence of these heat loss considerations.

When running the thermal simulation, to reduce required memory and computational time, the time points feature within Abaqus was used. Time points use a list of specified values to create a bespoke range of times within the analysis at which to output the solution. Note: time points are different to time increments, many increments may occur but only certain ones are written to the output file. It is used as a tool to reduce computational space required and does not affect the accuracy of the analysis. The total step time was taken as the entire build time and a single time point was used to output the thermal history for each elemental layer (e.g., each simulated, aggregation of actual build layers). A total of 79 frames were outputted within the thermal analysis.

Mechanical Simulation

Although not a numerical requirement, the mesh geometry is automatically copied from the thermal model and used in the mechanical model. This ensures mesh compatibility and minimizes any further interpolation required between nodes. The thermal elements are automatically replaced with 3D stress elements in the App. Similarly, the heat transfer step in the original heat transfer analysis was replaced with a single, static general step. The same step time and time points sequence were used, although a smaller time increment was used to ease convergence. Further details on this can be found in references [36, 37].

The thermal history was imported into the mechanical model as a predefined field. Temperatures are automatically mapped between the thermal analysis and structural analyses, thermal expansion then drives deformation. These values essentially act as the load within the analyses and can be used to determine the build-up of residual strain within the model.

The initial temperature of the build part represents a relaxation temperature (not room temperature) above which thermal straining induces negligible thermal stress. Upon material activation, it represents the temperature from which the initial thermal contraction occurs. In this analysis, the initial temperature of the bridge is set to 750 °C. The initial temperature for the build plate is 80 °C.

The thermal results for each increment during the previous transient heat transfer analysis are applied to the structural analysis as predefined fields. Abaqus automatically maps the nodal values of temperature by interpolation (both in space and time) of the previous results.

In addition to this, the base of the substrate is assigned an encastre boundary to ensure no movement. Therefore, we can determine from this that any displacement caused within the build is as a direct result of the build-up of residual stresses.

Results and Discussion

The thermal simulation, though done on a coarse scale, shows reasonable accumulation of heat among the layers and accounts for the cooling period between layers.

The results of the mechanical simulation are shown in Fig. 6. Figure 6a shows the contour plot of the residual strain in the x direction. Likewise, Fig. 6b shows the simulation result of the residual strain along the z direction.

Fig. 6
figure 6

Results of the mechanical simulation a EE11 and b EE33

When compared to the results produced by NIST [31], a few main features show strong agreement between the simulation and experimental results.

Firstly, in a comparison of the simulation and experimental contour plots in the x direction (Fig. 6a and Fig. 7a), one can see that both results show tensile strains across the main body of the part, with compressive strains being seen on the left hand side of each leg. Furthermore, there appears to be an accumulation of compressive strains around the peaks where the legs join. The magnitude in both cases is on a very similar scale, showing values in the region of 3.5e-03 and − 3.5e-03. Another prominent feature that can be seen in both contour plots is the change from tensile to compressive strains within the middle region of the far right support.

Fig. 7
figure 7

XRD measurements from NIST a EE11 and b EE33 [38]

The comparison of the contour plots in the z direction (Fig. 6b and Fig. 7b) again shows a very good level of agreement. Both contours show almost a reversal image to that in the x-direction. Compressive strains can be seen throughout the body of the build with tensile strains down the length of some of the legs, in particular the wider ones. Again, magnitudes are very similar, with a peak value of 3.652e-3 being shown in the simulations and 3.5e-3 in the experimental data.

Figure 8 and Fig. 9 show the comparison between X-ray diffraction measurements and simulation result of the strain in the z direction along the z = 10.75 mm and z = 2.75 mm paths in the center plane. The predicted results show good correlation with the benchmark test data.

Fig. 8
figure 8

XRD measurement and simulation result of strain z along the z = 10.75 mm path

Fig. 9
figure 9

XRD measurement and simulation result of strain z along the z = 2.75 mm path

It is suggested by the results above that the model implemented here has delivered a reasonably accurate prediction of the residual strains within the build. The magnitudes of the strains are correct across the plane examined and the distribution of the strains shows visual agreement between the two sets of contour plots. Additionally, Figs. 8 and 9 show very clearly an agreement in measurements throughout the build height of the part.

Conclusion

Within this study, a single cantilever was simulated and experimentally built and tested. This was undertaken as part of the Additive Manufacturing Benchmark 2018. Experimental data was supplied courtesy of the NIST laboratories.

A thermo-mechanical sequentially coupled analysis was undertaken in order to predict the residual strain values within the selectively laser melted part. The results show a very good level of agreement with the XRD measurements taken of the built part in this blind benchmark. From this, we can conclude that the simulation technique developed within this document is a sufficient method for the prediction of residual strain measurements. The model described in this paper received joint first prize in the residual elastic strain category of the NIST AM Benchmark 2018.