Thermophysical Phenomena in Metal Additive Manufacturing by Selective Laser Melting: Fundamentals, Modeling, Simulation and Experimentation
Among the many additive manufacturing (AM) processes for metallic materials, selective laser melting (SLM) is arguably the most versatile in terms of its potential to realize complex geometries along with tailored microstructure. However, the complexity of the SLM process, and the need for predictive relation of powder and process parameters to the part properties, demands further development of computational and experimental methods. This review addresses the fundamental physical phenomena of SLM, with a special emphasis on the associated thermal behavior. Simulation and experimental methods are discussed according to three primary categories. First, macroscopic approaches aim to answer questions at the component level and consider for example the determination of residual stresses or dimensional distortion effects prevalent in SLM. Second, mesoscopic approaches focus on the detection of defects such as excessive surface roughness, residual porosity or inclusions that occur at the mesoscopic length scale of individual powder particles. Third, microscopic approaches investigate the metallurgical microstructure evolution resulting from the high temperature gradients and extreme heating and cooling rates induced by the SLM process. Consideration of physical phenomena on all of these three length scales is mandatory to establish the understanding needed to realize high part quality in many applications, and to fully exploit the potential of SLM and related metal AM processes.
Additive manufacturing (AM) offers the opportunity to produce parts with high geometric complexity, without the requirement for dedicated tooling. AM processes for polymer parts are quite well established, yet melt-based processes for AM of metals still exhibit severe practical challenges, many of them resulting from the high melting temperatures of metals and their relatively low viscosities Lipson2013 (). In selective laser melting (SLM) of metals, the arguably most prominent representative of powder bed AM methods, a 3D manufacturing task is digitally segmented into thin 2D layers Gibson2010 (). A solid part is simply formed by selectively melting pre-defined contours in successive layers of powder using a focused laser beam. After one layer of powder has been scanned, the regions melted by the laser form the cross-section of the final part. Subsequently, the underlying build platform is lowered down and a further layer of powder is deposited by means of a powder coater mechanism. This procedure is successively repeated until the final 3D geometry is completed and the remaining unfused powder is then removed (see Figure 1). In this context, the so-called build direction denotes the direction normal to the powder bed. Further, the orientation of the part geometry with respect to this (vertical) build direction has crucial influence on the resulting part properties Hanzl2015 ().
SLM of metals offers significant advantages including near-net-shape production without the need for expensive molds or time consuming post-processing, high material
utilization rate and highest production flexibility Thijs2010 (); Stucker2011 (). Most essentially, the layer-wise production leads to nearly unlimited freedom of design, which enables the generation of highly complex geometries that cannot be obtained by conventional manufacturing processes. This paradigm shift in mechanical design allows to integrate complex substructures such as lattice-based geometries enabling lightweight yet sufficiently stiff components. Possible fields of application include aerospace or medical engineering and basically all industries requiring highly complex and individualized parts (see Figure 2).
The overall SLM process is highly complex and governed by a variety of (competing) physical mechanisms. The most important effects and physical phenomena occurring in the powder bed, the melt pool and the solidified phase of typical SLM systems are summarized in the following three subsections, as well as in Figure 1. A typical schematic parametrization of the powder layer and laser beam as well as a typical local scan pattern are shown in Figure 3.
1.1 Physical phenomena within the powder bed
The incident laser beam is a collimated, polarized, monochromatic electro-magnetic wave, with a wavelength in the range of . The spatial power density distribution of incident radiation on the powder bed is commonly assumed to follow a Gaussian distribution, with the associated value typically being taken as the laser beam spot size. Typical values for the nominal laser power and the employed laser beam velocities are in the range of and Thomas2016 (). The effective laser beam absorption within the powder bed is governed by multiple reflections of incident laser rays within the open-pore system of the powder bed, each with partial absorption of the incident radiation. The laser beam can penetrate to considerable depths, which can even reach the range of the powder layer thickness Gusarov2009 (); Gusarov2005 (). Thus the net absorptivity of powder beds is considerably higher than the value known for flat surfaces and, moreover, the laser beam energy source must to be thought of a volumetric heat source distributed over the powder bed thickness, as opposed to a surface heat source. The factors influencing overall absorption and local energy distribution are numerous, including the laser beam power, wavelength, polarization, angle of incidence, powder temperature, surface roughness, surface chemistry (e.g. oxidation) and contamination Boley2015 (); Khairallah2014 (). Issues of powder bed morphology, determined by particle shape, size distribution and packing density, are also central to radiative transfer. A further important factor is given by the intra- and inter-particle heat transfer within the powder bed. The inter-particle heat transfer is typically governed by the gas in the powder bed pores, with commonly negligible overall conductivity contributions from particle-to-particle contact points as long as loose, i.e. not mechanically compressed, powder layers are considered. Consequently, the thermal conductivity of loose powder is comparable to the conductivity of gas and by orders of magnitude smaller than the conductivity in the solidified phase Rombouts2005 (). Also when considering the intra-particle heat transfer, it can be observed that the time scales governing this process are typically larger than the time scales governing particle melting. In other words, under typical SLM process conditions, there is not enough time for conductive homogenization of non-uniform energy and temperature distributions across the powder bed but also across individual particles. As consequence, partially molten particles may cause defects such as pores or inclusions Boley2015 (). A further characteristic of the SLM process is that the powder bed, when considering it as a homogenized continuum, shows mesoscopic heterogeneities (in form of individual particles) that are in the same order of magnitude as relevant macroscopic process length scales such as powder layer thickness and laser beam spot size. Typically, powder particle sizes in the range of , layer thicknesses in the range of and laser beam spot sizes in the range of are employed. Also, in standard linear scan patterns the distance between two successive laser tracks, denoted as hatch spacing , is an important process parameter. This is typically chosen in the range of such that a sufficient overlap and remelting between two subsequent tracks is guaranteed (Figure 3). A good overview of the typically applied range of these process parameters is also given in Thomas2016 ().
The large size of individual powder grains as compared to powder layer thickness and laser beam spot size typically leads to non-uniform energy distributions, across the entire powder bed but also across individual particles, which may have considerable influence on the resulting melting behavior and melt pool hydrodynamics. Furthermore, these comparatively large heterogeneities cause differences in the resulting temperature fields and melt track shapes when considering different samples of stochastically equivalent powder layers (i.e. identical particle shape and size distribution, powder layer density and thickness etc.). Consequently, the variance of process results due to the stochastic nature of the powder layer is considerably greater than for comparable processes such as laser beam welding (LBW, see below). Besides the energy input by the laser beam, also possible energy losses in form of thermal radiation emission, thermal convection or heat conduction from the solidified material to the underlying built platform play an important role in the overall SLM process. For further information on powder bed radiation and heat transfer in the context of SLM, the interested reader is referred to Rombouts2005 (); Gusarov2005 (); Gusarov2008 (); Boley2015 ().
1.2 Physical phenomena within the melt pool
As soon as the melting temperature is reached at local positions on the powder grain surface, the phase transition from solid to liquid as well as the formation of a melt pool and, ideally, a continuous melt track is induced. Driven by surface tension and capillary forces tending to minimize surface energy, a coalescence of individual melt drops and a reshaping of the resulting melt pool is initiated Korner2011 (); Korner2013 (). In addition, the wetting behavior of the low-viscosity melt on the underlying substrate, i.e. the solidified material formed by previous layers, and surrounding powder grains influences the resulting melt pool shape, continuity and adhesion to the previous layer. The wetting behavior crucially depends on the material, temperature, surface roughness, and surface chemistry. Oxidation on the powder grain or substrate surfaces - either due to contaminated primary powder material or due to thermally induced oxidation during the process - is known to considerably decrease the wetting behavior of the melt which might result in instable, balled melt pools and rough surfaces, pores or delamination due to insufficient layer-to-layer adhesion Das2003 ().
The prevalent length and time scales essentially determine which physical effects govern the process and which are negligible. Typically, viscous and gravity forces can be considered as secondary effects while surface tension and capillary forces, wetting behavior but also inertia effects are the primary driving forces that influence the melt pool dynamics and shape as well as the surrounding powder morphology by attracting or rejecting individual grains Korner2013 (). The heat transfer within the melt pool is governed by convection rather than by heat conduction, with Marangoni convection, i.e. melt flow from hot to cool regions induced by temperature-dependent surface tension, playing a prominent role in this process Khairallah2016 (). Depending on the amount of absorbed energy density and the surrounding atmospheric pressure, the peak temperature within the melt pool might exceed the boiling temperature and considerable material evaporation may take place Gusarov2007 (). The evaporation itself as well as the gas flow induced by evaporation may influence the melt pool thermo-hydrodynamics and the overall process as consequence of an evaporative mass loss and additional cooling, of a recoil pressure considerably distorting the melt pool surface and representing a means of transport for potential pollutants, of melt drops spattered out of the pool and even of powder particles ejected away from the direct vicinity of the laser beam Khairallah2016 (). Driven by the so-called keyholing mechanism, cavitation resulting from evaporated material might considerably contribute to the overall material porosity or even burst the surrounding solidified material due to thermal expansion of trapped gas Liu2016 (). As soon as the melt pool has solidified, the evolution of possible defects on the mesoscopic scale is virtually established. A discussion of the physical phenomena governing the melt pool thermo-hydrodynamics can e.g. be found in Das2003 (); Kruth2005 (); Khairallah2016 (); Korner2013 (); Lee2015 ().
1.3 Physical phenomena within the solidified phase
With the solidification of the melt pool, the development of the metallurgical microstructure, crucially determining the macroscopic properties of the final part, begins. The evolution of the solid-phase microstructure characterized by grain size, grain shape (morphology) and grain orientation (texture) are governed by the prevalent spatial temperature gradients, the cooling rates, as well as the velocity of the solidification front Glcksman2010 (). In SLM processes, two regimes can be distinguished: The first regime is given by the temperature field in the direct vicinity of the laser beam, in the so-called heat affected zone (HAZ) Yadroitsava2015 (), which is controlled by highly complex mechanisms such as radiation absorption and heat conduction in the powder bed as well as convective heat transfer within the melt pool with all the individual physical phenomena and process parameters of influence as discussed in the two previous paragraphs. The material in this region is subject to a rapid heating above melting temperature due to the absorption of laser energy by powder grains, a high velocity of the melt pool front induced by the laser beam velocity as well as a rapid solidification of the molten material after the heat source has moved on, which is a direct consequence of the large ratio of solid material to hot molten material. These pronounced non-equilibrium conditions lead to meta-stable microstructures and compositions of the resulting phases, as well as smaller grain sizes - which typically result in higher material strengths - as compared to traditional melt-based manufacturing processes such as casting Herzog2016 (); Hanzl2015 (). A second regime of thermal evolution is prevalent in previously deposited material layers located below the current layer and further away from the heat source. Each solidified location experiences repeated heating and cooling cycles with decreasing amplitudes as the laser processes adjacent scan tracks, and as it processes consecutive new layers. The heat transfer in this regime is rather determined by global part properties, e.g. the global laser beam scanning strategy, the build direction, the fixation of the part on the built platform, the temperature of the built platform but also by the part porosity and the metallurgic microstructure distribution itself, which both influence the (effective) thermal conductivity.
Of course, microstructure also depends on the specific part geometry, e.g. due to heat flux concentration at the transition region from bulk material to slender columns or thin walls (see e.g. Figure 1), which are surrounded by the low conductivity unfused powder Mercelis2006 (); Lu2015 (). Also the evolution of columnar grain structures oriented in direction of the main temperature gradients, usually in build direction, is typical for SLM processes and often yields a strongly anisotropic macroscopic material behavior with higher material strength in the build direction Herzog2016 (); Gong2015 (). With increasing distance from the top powder layer, the maximal temperature values and gradients experienced by a material layer during the repeated thermal cycles decrease and, similar to a heat treatment, these cycles might lead to a coarsening of the microstructure, a reduction of brittle non-equilibrium phases and, consequently, to more ductile material characteristics. As consequence of these effects, the final part will typically exhibit a change of microstructure in built direction: On the one hand, the initial creation of fine grain structures and non-equilibrium phases will be more pronounced in the first material layers deposited in the direct vicinity of the build platform where higher thermal conductivity and faster cooling rates are prevalent. On the other hand, these initially deposited powder layers are exposed to the heat treatment of repeated heating and cooling cycles for longer times, which might lead to longer evolution times for solid phase transformations and grain coarsening.
Apart from the microstructure evolution considered so far, the high temperature gradients in the direct vicinity of the melt pool, but also at locations of heat flux concentration, are giving rise to considerable thermal strains induced by a successive thermal expansion and shrinkage of material. These thermal strains result in thermal stresses within the kinematically constrained SLM part. The magnitude of these stresses is essentially determined by the underlying solid microstructure and the resulting macroscopic material behavior. A ductile material can compensate these thermal strain variations by means of local plastic flow. On the contrary, brittle material behavior fostered by small grain sizes, the existence of certain non-equilibrium phases or a local segregation of alloying elements Thijs2010 () might result in cracks at locations of stress concentration such as residual pores or inclusions. The complexity of this thermo-mechanical coupling is further increased by the fact that the microstructure does not only influence the amount of residual stresses, but that also the prevalence of residual stresses, which are for example known to mechanically stabilize the metastable austenitic phase in steels Herzog2016 (), will influence the evolution of the microstructure. In order to reduce residual stresses on large surface areas, these surfaces are typically subdivided in smaller islands that are completed successively. Figure 3 (right) illustrates the processing of a single island consisting of contour scan and infill hatching. During the SLM process, the amplitude of residual stresses might also decrease since stress relaxation is likely to occur during the repeated heating and cooling cycles at lower temperature levels. Furthermore, annealing is typically applied to the final part before removing support structures in order to relief residual stresses. While a neglect of support structures would result in reduced residual stresses, this advantage has to be paid for by dimensional warping as often observed in parts made by SLM. For further information on material aspects and microstructure evolution in the context of SLM, the interested reader is referred to King2015 (); Bourell2016 (); Hebert2016 (); Herzog2016 (); Sames2016 (). Exemplary references focusing on the investigation of the residual stresses resulting during the SLM process are Hodge2014 (); Hodge2016 (); Denlinger2014 (); Riedlbauer2014 (); Yadroitsava2015 (); Lu2015 ().
All in all, it can be concluded that the entire thermal history between solidification and cooling down to the ambient temperature, governed by many heating and cooling cycles at different temperature levels and time scales, considerably determines the resulting metallurgic microstructure as well as the macroscopically observable material properties such as ductility, micro hardness, yield strength or tensile strength as well as their spatial distribution in a possibly inhomogeneous and anisotropic manner. On the one hand, a restriction of the overall process parameters is required in order to avoid undesirable material characteristics and possible part defects such as excessive residual stresses, dimensional warping, crack propagation or delamination of layers, effects which might destroy the SLM part during the build process or at least reduce the mechanical resilience of the final part considerably. In order to fully exploit the efficiency potential of SLM, this result has to be achieved without the need for time- and cost-intensive post-processing as required by alternative processes such as selective laser sintering (SLS, see below). On the other hand, the flexibility of the SLM process offers the unique opportunity to optimize process parameters in order to manufacture parts with prescribed inhomogeneous and anisotropic microstructures and macroscopic material properties in a controlled manner, contributing to the paradigm shift in design enabled by SLM.
1.4 Differentiation of related additive manufacturing processes for metals
References and comparisons to other powder-bed metal AM processes such as electron beam melting (EBM) and selective laser sintering (SLS) will be useful from time to time in the foregoing discussion Kruth2005 (); Kruth2007 (); Gong2014 (); Gong2014a (); Herzog2016 (). Similar to SLM, the EBM process represents a powder bed-based additive manufacturing process where pre-defined contours are selectively melted in successively deposited powder layers. While SLM applies a laser beam as energy source, EBM is based on an electron beam. EBM is only applicable to electrically conductive materials and has to be performed in a near-vacuum environment in order to avoid an interaction of the electron beam with surrounding gas molecules. On the other hand, SLM is also suitable for dielectric materials and has to be performed in an inert gas atmosphere in order to prevent surface oxidation of (metallic) powder and substrates. Also the way of energy transfer into the powder bed is different: When considering one individual powder grain, the laser beam radiation is (in good approximation) absorbed at the powder grain surface while in EBM, electrons penetrate the powder grain surface, a process in which kinetic energy of the electrons is dissipated leading to melting of the grain. However, the energy of the electron beam is typically completely deposited to the powder grain of first incidence. On the contrary, only a part of the laser beam energy is directly absorbed at the powder grain of first incidence while the remaining part is reflected, leading to considerably higher powder bed penetration depths due to multiple reflections in the open pore system provided by the powder bed. Eventually, the laser beam in SLM is suitable for smaller, more focused spot sizes as well as smaller powder grains sizes, which enables finer geometrical resolutions and an improved surface quality. The higher concentration of incident energy in SLM might yield higher local cooling rates and, thus, smaller grain sizes of the resulting metallurgical microstructure.
Similar to SLM, SLS is based on a laser beam energy source. However, in contrast to SLM and EBM, the particles within the powder bed are not fully molten. In this context, solid and liquid state sintering can be distinguished. Solid state sintering is governed by thermally activated diffusion of atoms at temperatures below the melting point leading to slowly growing necks between adjacent powder particles. However, the slow time scales governing this process and the resulting low output rates make it often infeasible from an economic point of view. On the contrary, liquid state sintering aims at a partial melting of powder. The molten liquid will typically spread almost instantaneously between the unmolten particles and acts as binder. In order to achieve defined ratios of molten and unmolten material, either a second material species with lower melt point is employed as binder, bi-modal powder mixtures of one and the same material are used such that the smaller grains melt earlier than the larger ones or the incident energy density has to be adjusted such that only the top surface of powder grains melts while the core remains in solid state. These sintering processes often result in so-called green parts with high porosity, which requires subsequent heat treatment processing.
Related processes such as directed energy deposition (DED) and wire-feed laser and electron beam additive manufacturing methods share some similarities with SLM and EBM, yet they generally operate on larger length scales ( scale molten trajectory). Also the processes of Laser Beam Welding (LBW) and Electron Beam Welding (EBW) will be considered at some points in this work, since some of the underlying physical mechanisms are similar to SLM and EBM. While SLM and EBM aim at additively generating solid material structures out of layers of loose powder, LBW and EBW intend the connection of individual solid parts by partially melting their contact surface.
1.5 Organization of this article
The remainder of this article is structured as follows: In Section 2, the governing physical mechanisms are further detailed, modeling as well as numerical (and analytical) solution procedures in the context of SLM processes are reviewed and findings derived by these approaches are discussed. Section 2.1 focuses on the modeling of radiation and heat transfer in powder beds. In Sections (2.2)-(2.4), a comprehensive overview of existing approaches classified by means of the three main categories of SLM modeling approaches found in the literature, namely macroscopic, mesoscopic and microscopic models, will be given. Section 3 focuses on experimental studies that are especially relevant to understand the thermophysical mechanisms. Following a similar structure, Section 3.1 deals with experimental approaches of powder bed characterization while Sections (3.2)-(3.4) refer to experimental investigations on effects visible on macroscopic, mesoscopic and microscopic level, thus representing the counterpart to the Sections (2.2)-(2.4) considering modeling and simulation. Finally, Section 4 provides recommendations concerning practical implementation of the SLM process and gives open questions and potentials for future process improvement.
2 Modeling and simulation approaches
Approaches to modeling of SLM can be classified in macroscopic, mesoscopic and microscopic models. Macroscopic simulation models typically treat the powder phase as a homogenized continuum resulting in efficient numerical tools capable of simulating the manufacturing of entire parts by SLM. Macroscopic models commonly aim at determining spatial distributions of temperature, residual stresses as well as dimensional warping within SLM parts. Mesoscopic models typically resolve individual powder grains and melt pool thermo-hydrodynamics in order to determine part properties such as adhesion between subsequent layers, surface quality as well as creation mechanisms of defects such as pores and inclusions. Microscopic models consider the evolution of the metallurgical microstructure involving resulting grain sizes, shapes and orientations as well as the creation of thermodynamically stable or unstable phases. The computational effort required for mesoscopic models currently limits the application of these models to single track simulations. However, the insights gained by mesoscopic models might serve as basis to further improve continuum models for powder and melt phase in macroscopic models. Similarly, microscopic models can readily address small areas of the part in the range of the powder layer thickness, and existing approaches are commonly limited to 2D. Nevertheless, they may serve as basis for the development of improved inhomogeneous and anisotropic continuum constitutive models whose quality is essential for the quality of simulation results derived by macroscopic simulation models. Macroscopic as well as mesoscopic models commonly require submodels for the radiation transfer into and the heat transfer within the powder phase. For this purpose, modeling approaches considering powder bed radiation and heat transfer are discussed in Section 2.1.
A challenge more severe than the actual derivation of physical models, is the application and development of powerful discretization techniques and numerical solution schemes in order to enable a robust and efficient computational solution, two key factors in the simulation-based characterization of SLM processes. While the finite element method (FEM), the finite difference method (FDM) as well as the finite volume method (FVM) represent spatial discretization schemes typically employed in the considered models, temporal discretization is almost exclusively based on explicit or implicit finite difference time integration schemes. For implicit schemes, the fully discretized problem is typically represented by a system of equations that is nonlinear in the unknown discrete primary variables, which requires the application of a nonlinear solver, e.g. a Newton-Raphson scheme. Depending on the characteristics of the nonlinear system of equations to be solved, convergence of the nonlinear solution scheme might not always be guaranteed. On the contrary, explicit schemes allow for a direct extrapolation of the known configuration at time to the unknown configuration at time . The resulting system of equations is linear in the discrete unknowns such that no iterative, nonlinear solution process is required. However, at least in the geometrically linear regime, implicit time integrators can be proven to be unconditionally stable, thus typically allowing for considerably larger time step sizes as compared to explicit schemes in order to preserve system stability, i.e. to keep the total system energy bounded during the simulation. Consequently, implicit schemes are favorable for problems that are dominated by a low frequency response, where the large time step sizes possible for these schemes are sufficient in order to resolve the low-frequent system answer. On the contrary, explicit schemes are rather suited to model high frequency responses and wave-like phenomena such as high velocity impacts. There, small time step sizes are required for both explicit and implicit schemes in order to accurately resolve the high-frequent system dynamics. Since the computational effort per time step is lower for explicit schemes and no nonlinear convergence has to be accounted for, these schemes are preferable in such scenarios. The stability-relevant time step sizes dictated for explicit schemes in the context of mesoscopic SLM models are typically considerably smaller than the time step sizes required to capture the relevant physical phenomena. Consequently, implicit schemes would have the potential to achieve substantial computational savings. However, the complexities arising from the multiple field and domain couplings as well as the geometrical characteristics prevalent in SLM make the implicit treatment of fully resolved mesoscopic models a challenging task.
A further important question to be asked from a numerical point of view is the way different physical fields (e.g. thermal and mechanical fields) and domains (e.g. powder, melt and solid phase) are coupled. Here, three concepts of treating phase boundaries shall be distinguished: The first category models a sharp boundary between the phases resulting in a jump of material properties and physical fields under consideration of displacement / velocity continuity conditions and mechanical equilibrium at the interface. Numerical realizations of such models are typically based on explicit interface tracking, e.g. via level set schemes, and discrete solution spaces allowing for discontinuity in the primary variables, e.g., the extended finite element method (XFEM). The second category of schemes still identifies the interface in an explicit manner, e.g., by means of an additional phase variable taking on the value for the one phase and for the other. However, the interface is not represented by a sharp boundary, but rather by a transition range of finite thickness characterized by phase variable values . An example for such a scheme is the volume of fluid method (VOF). A third category is given by methods that do not introduce additional variables in order to distinguish the phases. These schemes typically treat the phase boundary implicitly by means of specific values of existing primary fields, e.g., by defining the phase boundary between liquid and solid as the isothermal contour representing the melt temperature . The rapid change of physical properties at the phase boundary is considered by these schemes in terms of high gradients in the applied material parameters. Thus, again, the interface is extended to a finite thickness. A too small transition region might lead to excessive gradients in the material parameters and, as consequence, to ill-conditioned discrete problems that are challenging to be solved numerically. The SLM modeling approaches reviewed in this work almost exclusively rely on the second and third category of approaches.
Numerical schemes can also be distinguished concerning the succession of solving different physical fields and domains. Namely, monolithic schemes solve the entire multi-physics problem at once, and partitioned schemes solve the different physical fields/domains subsequently. The partitioned schemes can be further subdivided into so-called iterative or strongly coupled partitioned schemes that iterate between the different fields several times within a time step until achieving the solution of the monolithic problem statement and so-called staggered or weakly coupled partitioned schemes that solve for each field only once per time step without additional iterations, eventually leading to a solution that differs from the monolithic one. Monolithic and strongly coupled partitioned schemes are typically combined with implicit time integrators because these coupling schemes can preserve the desirable stability properties of implicit integrators at large time step sizes. On the other hand, weakly coupled schemes are often combined with explicit time integrators since these schemes also lead to low computational costs per time step. However, for weakly coupled partitioned schemes, even when combined with implicit time integrators on the individual fields, the stability requirement often becomes very restrictive typically dictating very small time step sizes and in some cases even prohibiting stability at all (see e.g. Forster2007 (); Causin2005 () for a discussion in the context of fluid structure interaction (FSI) simulations). Since time step size correlates with computational costs, these considerations are of high practical interest.
2.1 Modeling of optical and thermal properties of powder phase
The optical and thermal properties of the powder bed, in combination with the laser characteristics, crucially determine the heat distribution in the powder bed and the subsequent melt pool dynamics. This section addresses approaches to modeling the optical and thermal properties of the powder bed as well as the radiation-dominated energy transfer from the laser beam source into the powder bed and the conduction-dominated heat transfer within the powder bed. The considered approaches are typically employed as powder bed submodels within the three main categories of macroscopic, mesoscopic and microscopic models. Often, the laser beam can be described in good approximation by means of a Gaussian power distribution Dong2009 (); Verhaeghe2009 (); Roberts2009 (). In a first step, the question of how the laser energy transfers into the SLM powder bed, a process that can be described via the principles of thermal radiation, geometrical optics and electro-magnetic wave propagation, will be considered. When it is not mechanically compressed after spreading, the powder bed has as high porosity as freely poured powder, which is in the range of for typical SLM powders Rombouts2005 (). As shown in Gusarov2009 (); Gusarov2005 (), laser radiation penetrates into powder through pores to a depth of several particle diameters because of multiple reflections Wang2002 (). This depth is comparable with the powder layer thickness. Thus, laser energy is deposited not on the surface but in the bulk of the powder layer. The incident laser energy consists of portions that are reflected, absorbed and transmitted when impinging upon the surface of powder particles. The emission of radiation from the powder particles themselves can often be neglected. In terms of radiation reflection, it can be distinguished between specular (mirror-like) and diffuse reflection, with the latter exhibiting a reflection intensity that is equally distributed over all possible directions. Concerning radiation transmission, typically the classification of opaque, semi-transparent or transparent particles is made.
The multiple reflection/absorption processes of the light in the powder bed is additionally influenced by powder characteristics such as mixture ratio, mean particle size and shape, size distribution, packing density, powder bed depth but also by the laser beam spot size. Also the polarization of the laser beam relative to the powder particle surfaces can lead to non-uniform energy absorption even across individual powder particles. Since the thermal conduction in the powder bed is typically weak due to the prevalent porosity and melting time scales are typically smaller than heat conduction time scales for individual grains, non-uniform energy absorptions across the powder bed but also across individual powder particles will in general have considerable influence on the melt pool shape and the properties of the solidified track. A further factor of influence is given by possible surface oxidation and contamination effects.
An extensive review on the topic of radiation transfer within heterogeneous/porous/dispersed media can be found in Viskanta2016 () as well as Baillis2000 (). The following section will mainly focus on approaches in the context of SLM. According to Baillis2000 (), the modeling approaches for radiation transfer in heterogeneous media can be classified as models based on a continuum formulation of the radiation transfer equation (RTE) well-known for homogeneous (participating) optical media, and models based on a discrete formulation of the RTE, which typically leads to ray tracing schemes. Both approaches are based on the simplifying assumptions of geometrical optics which considers light rays that: propagate in rectilinear paths as they travel in a homogeneous medium; bend, and in particular circumstances may split in two at the interface between two dissimilar media; follow curved paths in a medium in which the refractive index changes; and may be absorbed or reflected.
2.1.1 Continuum model for powder bed radiation transfer
The general radiation transfer equation (RTE) known for homogeneous media and underlying the homogenized, continuum models for powder bed radiation transfer is (see also Baillis2000 (); Gusarov2005 (); Chandrasekhar2013 ()):
The RTE describes the rate of the direction- and position-dependent radiation energy flux density based on an energy balance considering radiation absorption (first term on the right-hand side), radiation emission (second term on the right-hand side) and scattering (first and third term on the right-hand side). In this context, represents a directional unit vector, is an infinitesimal solid-angle increment whose magnitude is per definition identical to an infinitesimal surface element on a unit sphere, and the product represents the vector of energy flux density transferred by photons in direction of the unit vector within a solid angle increment at position . The constants and are the scattering and absorption coefficients, which are, however, often replaced by the alternative constants extinction coefficient and albedo (portion of reflected light) . The mapping is the (normalized) scattering phase function stating the probability that radiation in -direction is scattered to the -direction and is Planck’s blackbody function representing radiation emission at a certain position:
Here, is the powder temperature, is the temperature of the ambient gas atmosphere and is the Stefan-Boltzmann constant. While the radiation transfer equation has originally been derived for homogeneous media, it is well-established to apply (1) also to heterogeneous systems, where the continuous quantities and parameters introduced so far have to be replaced by their effective counterparts, determined via spatial averaging over a reference volume that has to be much greater than the length scales of prevalent heterogeneities (e.g., of individual particles in powder beds). When considering radiation heat transfer in powder beds, the extinction coefficient is typically determined by the structure of powder, i.e., by size and shape of particles and by their arrangement, but is independent of optical properties of the material of particles like reflectance. On the contrary, the scattering characteristics such as the albedo and phase function commonly depend on reflective properties of the material of particles.
In the following, a coordinate system is chosen such that the positive -axis points into powder layer thickness direction with representing the upper surface of the powder bed and the boundary between powder bed (with thickness ) and the underlying solidified phase (substrate). The net radiation heat flux density , following the typical definition as heat flux per unit area measured in , results from the energy flux density per unit angle increment via integration over all directions of incident radiation. Expressed by means of solid-angle increments this is equivalent to an integration over the surface of a unit sphere. Assuming that the and components of the heat flux density resulting from a solution of (1) are negligible, also the incident power density per unit volume can be determined based on the following relations:
Gusarov et al. Gusarov2005 () considered the general problem of normal incidence of collimated radiation on a thin powder layer consisting of opaque (i.e. optically non-transparent) particles of arbitrary shapes and dimensions placed on a reflecting substrate in order to derive resulting net absorptivities and deposited energy profiles. It is assumed that the particle size is sufficiently small as compared to the laser spot size and powder bed thickness as to allow for the employed homogenization procedure. Further, in SLM the particle size is typically much greater than the wavelength of radiation, thereby justifying the applicability of the geometrical optics theory. The contribution representing the emission from individual powder grains has been neglected due to the comparatively high energy density of the incident laser radiation. Gusarov et al. Gusarov2005 () propose a statistical model based on the powder porosity and specific powder surface areas in order to determine the model parameters required in (1). This statistical model accounts for dense powder beds where particles touch each other and radiation transfer cannot be treated as scattering by independent particles, an effect that is referred to as dependent scattering (see e.g. Drolen1987 (); Chen1963 (); Kamiuto1991 (); Jones1996 (); Kamiuto1990 (); Tancrez2004 ()). Moreover, an analytical solution for the radiation transfer equation has been derived via the so-called two-flux method with the boundary conditions of normally collimated incident flux and specular reflection at . Among others, the additional assumption of the laser spot size being much larger than the absorption depth has been made. Thus, the radiation transfer problem could be considered as 1D-problem in direction with an incident flux .
In Figure 4, the resulting energy deposition profiles derived by the RTE solution as proposed in Gusarov2005 () and by means of ray tracing simulation are compared for three examples considering either WC-Co or Fe-Cu mixed powders with particles sizes (Fe and WC), (Cu) and (Co) as well as laser beam wave lengths of or . Accordingly, these two modeling approaches are in good agreement since a sufficiently deep powder bed and a sufficiently large laser beam size has been chosen. The total/integrated laser energy absorbed in a thin powder layer on a reflective substrate increases with its thickness while the deposited energy density per unit volume decreases. The absorption depths in Figure 4 seem to be comparatively high as compared to typical values known from SLM powder layers, which might be an indication for rather low packing densities considered in this study. Subsequently, the model of Gusarov2005 () has e.g. been applied in the independent works Verhaeghe2009 (), Khairallah2014 () and Hodge2014 (). In Gusarov2008 (), the model of Gusarov2005 () has been further extended and refined. There, it has been shown that the application of the model derived in Gusarov2005 () is only reasonable if either one of the two phases (powder particles and pores) are opaque, one phase has a much larger volume fraction, or the radiation properties of the two phases are at least similar. On the contrary, an application of this model, which simply averages the radiation intensities across both phases, to problems where the two phases capture comparable volume fractions but show considerably different optical properties at the same time (e.g., mixed powders), seems unjustified.
For mixed materials, a model based on partial homogenized radiation intensities and denoted as Vector Radiation Transfer Equation (VRTE) has been proposed in Gusarov2008 (). In this model, the partial values are obtained by averaging over each individual phase. It consists of two transport equations for the partial homogenized radiation intensities. The equations are similar to the conventional RTE, but contain additional terms that take into account the exchange of radiation between the individual phases. The structure of the medium is specified by volume fractions of phases and by specific surfaces of phase boundaries. The optical properties are described by refractive indices and absorption coefficients of phases. It has been verified that the vector RTE model reduces to the conventional RTE model in case one of the two phases is opaque or one phase prevails in volume. Compared to the RTE in Gusarov2005 (), an analytical solution was no longer achievable for the VRTE Gusarov2008 (). Instead, a discrete-ordinate method Carlson1965 () was exploited.
2.1.2 Ray tracing model for powder bed radiation transfer
As pointed out in Boley et al. Boley2015 (), the assumption of a RTE continuum model Gusarov2005 () is questionable for very thin, low-porosity metal powder layers with a layer thickness in the range of a few powder particles and/or a laser spot size comparable to the size of the powder particles. In such scenarios, which are often prevalent in SLM processes, the averaging process underlying the RTE continuum models lead to considerable model errors as compared to the exact solutions, which typically yield an energy absorption that is highly non-uniform with respect to the spatial position of the laser beam. Ray tracing modeling is one possible approach to resolve such heterogeneities. However, it shall be emphasized that depending on the spatial resolution required by the physical question to be answered (e.g. when global residual stress distributions but not the specific locations of individual pores are of interest), also the application of RTE continuum models might be reasonable, which typically require a considerably lower computational effort as compared to ray tracing simulations.
In ray tracing simulations, the total energy emitted by the laser beam in a certain time interval is represented by a discrete ensemble of rays with defined spatial position, orientation and energy. The position and energy associated with the individual rays is typically chosen such that the overall energy emission but also the spatial energy distribution resulting from the entire ensemble equals the corresponding characteristics of the laser beam (e.g. a Gaussian energy distribution). After defining the individual rays, the path of each ray is traced until striking an obstacle (powder particle). Based on the optical properties of the obstacle surface, part of the ray energy is absorbed whereas the remaining part of the ray energy is represented by a reflected ray with defined energy and orientation, which will further be traced trough the powder bed (Figure 5(a)). Commonly, several reflections are considered for each ray until the remaining energy drops below some predefined threshold. The individual ray energy contributions absorbed by each particle are accumulated during the simulation. Also the ray-tracing model is based on the principles of geometrical optics, thus, the corresponding requirements mentioned above have to be fulfilled. In particular, the particle radius has to be considerably larger than the laser wavelength.
Ray tracing is a well-established approach for studying radiation problems in countless scientific disciplines. Several contributions have applied this approach in order to study the radiation transfer trough porous media and packed beds Chan1974 (); Yang1983 (); Singh1991 (); Singh1992 (); Singh1994 (); Argento1996 (). In one of the first ray-tracing approaches in the context of SLM/SLS, Wang et al. Wang2000 (); Wang2002 () studied a powder mixture consisting of two species of spherical particles differing in particle size and material. In Wang2002 (), the dimensions of the sintering/melting zone are estimated by determining the accumulated energy absorbed by each particle. However, no inter- or intra-particle mechanisms of heat transfer (such as heat conduction or convection) besides radiation has been considered. Based on a powder bed consisting of Tungsten Carbide (WC) particles and Cobalt (Co) particles and typical SLS system parameters, the ray tracing simulations as well as accompanying experiments came to the consistent result that of the absorbed energy is concentrated within a depth of from the powder bed surface. The resulting spatial distribution of the absorbed energy , normalized by the total incident energy , is plotted over the powder bed depth coordinate in Figure 6. While the highest amount of primary laser beam radiation obviously arrives directly at the top surface of the powder bed, the maximum of the absorbed energy distribution occurs slightly beneath the top surface since at this location also secondary radiation contributions, stemming from multiple reflections in the powder bed, are prevalent.
In Boley2015 (), the model of Wang2000 (); Wang2002 () has been further refined by additionally considering the angular- and polarization-dependence of the laser absorption. Thereto, the absorptivity of the laser beam components in perpendicular (index ) and parallel (index ) polarization states have been distinguished according to the Fresnel formulas (see Landau1984 ()):
This description takes into account the model of laser beam radiation as an electromagnetic wave with speed and direction of propagation . This wave is assumed to be (linearly) polarized with the constant unit vector representing the direction of the electric field normal to . In (4), the complex refraction index , describing the speed of electromagnetic wave propagation (real part ), absorption of the electromagnetic wave (imaginary part ) and the angle of incidence between and the normal vector onto the powder particle surface have been employed. The perpendicular (index ) and parallel (index ) components of the polarized electromagnetic wave are typically determined by means of a projection of the polarization vector into a reference plane spanned by the vectors and . In Figure 5(b), the resulting spatial absorptivity distribution across one isolated sphere of stainless steel irradiated by a horizontally polarized laser beam is depicted. Correspondingly, the distribution of energy absorption can vary considerably across one individual particle.
The high importance and practical relevance of these fluctuations shall be illustrated by the following time scale estimate. The time scales governing the homogenization of energy on a sphere with radius R due to thermal conduction can be estimated by , where is the thermal diffusivity of the metal. The time scales governing the melting of a spherical particle can be approximated by , where is the melting enthalpy per volume, is the laser irradiance, and is the flat-surface absorptivity. For parameter values typical for SLM processes, the melting time is often considerably shorter as compared to the diffusion time. Consequently, non-uniformity of energy absorption results in only partial melting of the particle, which might in turn considerably influence melt pool dynamics or creation mechanisms of pores Boley2015 (). However, in cases where melting takes place on larger time scales than heat conduction, the initial energy and temperature distributions might homogenize before melting occurs, and, consequently, approaches considering a constant energy distribution across a particle (see e.g. Wang2002 ()) might yield reasonable results. The fact that thermal conductivity of loose powder is typically governed by the gas filling the powder bed pores Rombouts2005 () yields even larger time scales for the inter-particle heat conduction.
According to Boley2015 (), there are two general factors of influence leading to absorption non-uniformity, one related to the non-uniformity of absorption within a single particle as considered above, and the other related to the non-uniformity of the powder bed. In order to investigate this second factor of influence, ray tracing simulations on laser tracks across powder beds (thickness ) resting on a solid substrate and consisting of two different powder mixtures, a Gaussian distribution (average size , Figure 7(a), top) and a bimodal distribution with maximal packing density (particle ratio :, Figure 7(a), bottom), have been conducted in Boley2015 (). For both powder mixtures, typical distributions of the total absorbed energy along the scan track length as illustrated in Figure 7(b) can be observed. While large laser beam sizes (see curve in Figure 7(b)) lead to a comparatively homogeneous spatial energy absorption, laser beam sizes in the range of the particle diameter (see curve in Figure 7(b)), lead to strongly fluctuating spatial absorption values. Furthermore, the overall absorptivity turned out to be higher for the bimodal distribution, although the degree of fluctuations was comparable for both mixtures. Again, these fluctuations are not considered in homogenized continuum models such as Gusarov2005 (). Thus, their applicability requires sufficiently small particle sizes as compared to laser spot size and powder layer thickness.
In Boley2016 (), the results derived in Boley2015 () have been extended to a considerable range of different materials and verified experimentally. Concretely, the simulated power bed absorptivity and the flat surface absorptivity of the considered materials have been compared. Plotting the powder bed absorptivity over the flat surface absorptivity has revealed a smooth functional relation between these two quantities as already observed in Gusarov2006 (); Gusarov2010 (). This interesting observation allows to extract the assumed powder bed absorptivity for a material with given flat surface absorptivity, a procedure that can be considered to be of highest practical relevance (Figure 8). In Zhou2009 (), a comparable ray tracing simulation model has been proposed in the context of SLS.
2.1.3 Continuum models for powder bed heat conduction
The macroscopic simulation models discussed later in Section 2.2 regard the powder bed as a homogenized continuum. Consequently, these approaches typically rely on a model for an effective, homogenized powder bed conductivity. On the contrary, mesoscopic simulation models as considered in Section 2.3 resolve individual powder grains by spatial discretization. Consequently, each powder particle can be considered as a separate solid body, which simplifies the powder bed heat conduction problem to the two subproblems of intra-solid heat conduction as well as inter-solid heat conduction across the contact interfaces.
Study of thermal transport in two phase systems has long been a topic of scientific inquiry, with Maxwell and Rayleigh each publishing models of the problem in 1873 (see Maxwell1873 ()) and 1892 (see Rayleigh1892 ()), respectively. These early contributions have been based on comparatively strong model assumptions. For example, Rayleigh’s treatment idealizes the geometry to regularly spaced spherical particles in non-contact within the dispersed phase. An extension of the spherical particles to a larger spectrum of geometric primitives has been proposed in the contribution by Bruggeman Bruggeman1935 (). In general, the aim of more recent work has been to further relax such assumptions and match the model to expected behavior in limit conditions, often in the context of interpreting experimental results Tsotsas1987 (). For example, an often cited model has been published by Meredith and Toblas to explain the conductivity of water-propylenecarbonate emulsions. This approach is assumed to yield more realistic predictions than Maxwell’s or Bruggeman’s models at high volume fractions of the dispersed phase Meredith1961 (). In the review by Tostsas and Martin Tsotsas1987 (), models are further divided by how they treat primary parameters of bulk conductivities and porosity, and secondary parameters including radiative transfer, pressure dependence, contact between particles and deformation as well as convective effects.
The effective thermal conductivity of loose metallic powders, i.e. of material with negligible particle-to-particle contact surfaces, is typically controlled by the gas in the pores Rombouts2005 (); Gusarov2007 (). On the contrary, in Gusarov2003 (), a model has been derived for the heat conduction in power beds in case of small but finite contact areas between the particles. Such a scenario typically arises during the sintering process in SLS, at the beginning of the powder particle melting process in SLM or already in the non-melted solid state of the powder bed in case it is not loose but pre-compressed.
Now, the effective thermal conductivity within a powder bed in direction of a given unit vector can be defined on the basis of the heat flux density and the temperature gradient according to:
For illustration, the effective thermal conductivity shall briefly be presented for two well-known, but conceptionally different, heat transfer models for heterogeneous / particle-based systems, the Maxwell model mentioned above (index M) and the Reimann-Weber model (index RW). According to these models, the effective thermal conductivities are:
The well-known Maxwell model describes the effective thermal conductivity of systems consisting of randomly packed spheres with conductivity embedded in a continuous medium with conductivity and is applicable for sufficiently small volume fractions of spheres. It can easily be verified that this model yields the expected results for limit conditions such as , or . On the contrary, the Reimann-Weber model considers the effective thermal conductivity between two isolated spheres (bulk material conductivity , radius ) exhibiting a circular contact area of radius . The packing densities prevalent in powder bed AM make the application of the Maxwell model to these processes virtually impossible. For that reason, Gusarov et al. Gusarov2003 () derived formulations for the thermal conductivity of regular and random packings of spherical particles based on models of the ”isolated-spheres” type comparable to the Reimann-Weber model. As such, the effective conductivity of ordered powder beds has been determined by summing up the conductivity contributions of all individual sphere-to-sphere pairs. In the case of random packings, the models of ”isolated-spheres” type have been homogenized by spatial integration and the geometrical quantities and are replaced by effective quantities such as the particle volume fraction , the mean coordination number (describing the average number of closest neighbors) and the average dimensionless contact size ratio leading to the following effective thermal conductivity (index G for Gusarov):
While a strong dependence on the contact size , at least qualitatively, has already been known from experiments, experimental characterization approaches typically consider only the powder density , but not the coordination number . However, the derivations in Gusarov2003 () recommended for future experiments to treat these two parameters independently since for randomly packed powder beds different coordination numbers, and according to (7) considerably differing effective conductivities, are possible for the same powder bed density . Besides modeling approaches for the determination of effective powder bed conductivities, often also experimentally determined effective conductivities are employed in macroscopic SLM models as discussed in Section 2.2. An extensive comparison of experimental and modeling approaches for effective powder bed conductivities can also be found in Yagi et al. Yagi1957 ().
In summary, ray tracing models are based on a comparatively simple theory, yield a high degree of detailedness in the obtained solutions, but also require considerable computational resources to be applied to scenarios of practically relevant size. On the other hand, heat transfer continuum models are based on a more complex theory, are in general not capable of accurately resolving details on the length scales of prevalent heterogeneities (i.e., on the length scale of individual powder particles in the case of SLM), but are considerably less expensive for computation. Some simple cases even allow for the derivation of analytical solutions, for example by Gusarov Gusarov2005 (). Importantly, ray tracing models require specification of the powder bed structure by particle shape, dimension, and their coordinates, which is a non-trivial problem itself and often requires additional assumptions. In principle, the distinction of these two model classes is comparable to the distinction of atomistic and continuum models in other physical disciplines such as electricity or mechanics. Unfortunately, in the characterization of SLM processes, one is often interested in length scales comparable to the powder layer thickness. Because powder bed heterogeneities in form of individual particles are on the same length scale, the application of continuum models can often only be considered as an estimate of the underlying radiation transfer processes. Moreover, the assumption of a laser beam penetrating into a powder bed by means of multiple reflection actually only applies for sintering or for SLM processes based on point-wise modulated scanning strategies as long as no melt pool is prevalent. When SLM is performed with linear/continuous scanning, however, the resulting melt pool front typically proceeds the laser beam position and, thus, the laser beam impacts directly on the melt pool and not on an undisturbed powder bed. In order to accurately resolve the impact of the laser beam on the melt pool, whose actual surface contour is typically considerably distorted and part of the solution itself, ray tracing techniques seem to be indispensable. Similar conclusions can be drawn by comparing mesoscopic and homogenized macroscopic models for the powder bed heat conduction. Because the local energy and temperature distribution within the powder bed but also within individual powder grains crucially influence the melting behavior on the mesoscopic scale and the creation of defects such as pores or inclusions, also sufficiently resolved mesoscopic models for the laser beam radiation transfer (e.g. via ray tracing) and the inter- and intra-particle heat conduction (e.g. by resolving these particles in the spatial discretization of the thermal problem) are mandatory if questions on this scale are relevant. However, when applying macroscopic SLM models, i.e. continuum approaches suffering from a powder homogenization error anyways, heat transfer continuum models appear to be a reasonable choice.
2.2 Macroscopic simulation models
Macroscopic simulation models in the context of SLM processes typically treat the powder phase as a homogenized continuum described by means of effective, i.e., spatially averaged, thermal and mechanical properties, without resolving individual powder grains. This homogenization procedure yields efficient numerical tools capable of simulating entire SLM builds of practically relevant size across practically relevant time scales. The thermal problem lies typically in the focus of interest when these models are applied. In some works, an additional solid-mechanics problem statement is considered aiming at the assessment of global residual stress distributions or dimensional warping effects based on the full thermo-mechanical interactions.
The thermal problem is commonly given by the following set of balance equations, boundary and initial conditions:
Here, (8a) represents a variant of the energy equation as typically employed in the macroscopic and mesoscopic SLM models on the problem domain within the considered time interval , where is the density, is the specific heat at constant pressure, is the temperature and represents the velocity field. In (8a), is the thermal conductivity tensor and represents a heat source term. The two terms on the left-hand side of (8a) represent the material time derivative of the thermal energy density constituted of the local and the convective time derivative. It has to be mentioned that the thermodynamically consistent statement of thermo-mechanical problems involving compressible materials would require additional thermo-mechanical coupling terms in (8a). However, most of the models discussed in the following typically assume either exact or approximate incompressibility and resign these additional terms. Typically, isotropic conductivity is assumed, which yields . Often, the source term is modeled based on a solution of the RTE according to (3), e.g., the one proposed in Gusarov2005 (). Equation (8b) accounts for essential boundary conditions with prescribed temperature on the boundary whereas (8c) represents natural boundary conditions with prescribed heat fluxes on the boundary typically accounting for thermal convection and radiation emission on the powder surface with normal :
Here, denotes the temperature of the ambient gas atmosphere, is the convection coefficient and the emissivity. The Stefan-Neumann equation (8d) describes the phase change due to melting at the interface between the powder phase and the melt pool. In this equation, the indices and refer to the temperature gradients and thermal conductivities in the solid and liquid phase, respectively. Moreover, is the melting temperature, is the latent heat of melting and is the normal vector of the solid-liquid interface defining the melt process as a pure interface phenomenon. Eventually, (8e) prescribes the initial temperature .
If the solid mechanics problem is also solved as well, the following system has to be considered in addition to (8):
Equation (10a) represents the mechanical equilibrium of linear momentum. In this equation, the Cauchy stress tensor is related to the primary displacement field and temperature field by constitutive parameters with contributions from elastic , plastic and thermal strains. The mechanical equilibrium of angular momentum is fulfilled by definition due to the symmetry of the Chauchy stress tensor. The total time derivative represents the material acceleration vector and is the vector of volume forces acting on the physical domain . Similar to the thermal problem, essential and natural boundary conditions are given by equations (10b) and (10c) prescribing displacements and tractions . Equation (10d) represents the requirement of displacement continuity and mechanical equilibrium at the relevant interfaces, e.g. the interface powder-melt or the interface melt-solid, characterized by a normal vector field . The superscripts and denote quantities on the two different sides of the interface. However, most of the considered macroscopic models do typically not resolve these interfaces in the sense of a sharp 2D interface with discontinuous material parameters, but are rather based on a smooth, homogenized transition between the different phases, with the latter variant being easier to realize numerically. Finally, the initial position and velocity field is given by equations (10e) and (10f).
The spatial discretization of the heat equation (8a) and the momentum equation (10a) is typically based on the finite element method (FEM), which requires transfer of these equations into the (equivalent) weak form. For time integration, explicit as well as implicit approaches can be found. The coupling between the thermal and the solid-mechanical problem is often realized in a staggered partitioned manner. In the following, some of the most important representatives of macroscopic models available in the literature will briefly be discussed. As long as macroscopic models are considered, melt pool dynamics are not explicitly resolved and porosity-dependent, effective material properties are employed. The principal strategies of deriving macroscopic models are comparable for SLS, SLM and EBM processes in many aspects. Thus, scientific contributions related to all of these processes will be considered.
In recent years, considerable research effort in the development of macroscopic, mesoscopic and microscopic models has been conducted at the Lawrence Livermore National Laboratory (LLNL), including a macroscopic homogenized continuum model for SLM processes described by Hodge et al. Hodge2014 (). In this model, the laser beam heat source term has been taken from Gusarov2005 (). Heat losses occurring for example due to thermal emission, evaporation or mass ejection are taken into account by reducing the nominal total laser power to an effective total laser power, . In order to describe the phase transition, two spatial phase function fields and with are introduced representing the powder and the consolidated phase., i.e. for pure powder and for pure consolidated material. Importantly, the consolidated phase considered in this model is assumed to have a vanishing porosity and captures both, the melting as well as the solidified phase. In other words, only the phase transition from powder to liquid during melting is explicitly considered by means of a phase function, while the phase boundary between melt pool and solidified material is accounted for by means of high gradients in the material properties at these spatial locations. All of the relevant thermal material parameters (e.g. conductivity, heat capacity etc.) are considered to be a function of both temperature and phase (thus, of porosity). In the spatial discretization process based on the finite element method (FEM) and the subsequent numerical implementation, the phase boundary between powder and melt pool is not resolved in a sharp manner. Instead, there is typically a small band of finite elements where both phases are prevalent, a fact that is reflected by non-vanishing values . Correspondingly, the Stefan-Neumann equation is not evaluated in its original form (8d) at a 2D interface, but rather in an equivalent form within 3D volume elements.
While the phase boundary between powder and melt pool is tracked by means of phase variables, the Stefan-Neumann equation (8d) (or a spatially average version of it) is not considered for the phase transition from melting to solid. Instead, a high value of the heat capacity is chosen in the range of the melting temperature in order to account for the released latent heat during the solidification process. Thus, the typical temperature plateau observable in enthalpy-temperature diagrams of elementally pure materials at the melting point is replaced by an interval of very slowly increasing temperature. From a physical point of view, this might rather be comparable to the melting behavior of alloys between liquidus and solidus temperatures and . From a numerical point of view, it can be interpreted as a (commonly employed) regularization of the constraint equation in (8d) of elementally pure materials in order to simplify the numerical solution process.
Hodge et al. Hodge2014 () considered the thermo-mechanical problem by coupling the thermal and mechanical simulations in an iteratively partitioned manner. There, an elasto-plastic material law, additionally accounting for thermal expansion and consolidation shrinkage, has been employed while the boundaries between powder, melt and solidified phase have been considered implicitly via strong gradients of the associated mechanical material parameters with respect to temperature and porosity. This means, the interfaces given in equation (10d) are not explicitly resolved in the mechanical simulation. Furthermore, this elasto-plastic material model has been applied to all three phases, i.e. to the powder phase, to the solidified phase, and also to the liquid phase. Thus, a detailed modeling of the melt pool fluid dynamics, as prevalent in the mesoscopic models discussed in Section 2.3, is missing. This statement applies in a similar manner to most of the macroscopic SLM models considered in this section.
The thermal simulations in Hodge2014 () confirmed the results of Gusarov et al. Gusarov2007 () for 316L stainless steel concerning maximal temperature, top surface shape as well as cross-section shape of the melt pool. It was observed that the melt pool shape gets narrower and longer with increasing scan speed and that at some point a concave region at the tail of the melt pool boundary occurs, an effect that can be explained by the higher thermal conductivity of the solidified material behind the melt pool as compared to the powder material to the sides (Figure 9). Nevertheless, when considering the similarities between the two models proposed in Gusarov2007 () and in Hodge2014 (), it has to be emphasized that both approaches are based on the same radiation transfer model for the incident laser energy, and, that the chosen submodel for the heat source might considerably influence the overall results since the thermal conductivity within the powder phase is strongly limited by the low conductivity of the atmospheric gas filling the powder pores.
Because the model in Hodge2014 () explicitly considered the powder-melting phase change via and , the inertia of the melting process could be visualized by means of a melt pool boundary lagging behind the isothermal of the actual melting temperature in the range of high scan velocities. Furthermore, Hodge2014 () investigated the melt pool behavior when scanning across overhangs supported by loose (poorly conductive) powder (Figure 10). Experimental results Craeghs2010 (); Wang2013 () could be confirmed that typically larger melt pool sizes, excessive overheating and evaporation occur when applying the same scan parameters in domains of low thermal conductivity such as overhangs or thin-walled features. Commonly, this problem is addressed by adapting the laser parameters and/or increasing the net thermal conductivity in this regions by means of additional support columns.
Further validation of the Hodge simulation via experiments is given in the follow-up work Hodge2016 (), which also investigates the effects of build orientation, infill scan strategy, and preheating on residual stresses and dimensional warping. Digital Image Correlation (DIC) was used to measure the displacements of 316L triangular prisms fabricated by SLM as test artifacts. Moreover, neutron diffraction experiments were used to measure the stress distribution. Figure 11 illustrates typical results for a triangular prism printed vertically, with the shortest edge of the triangle forming one edge of the rectangular contact area with the build platform. The prisms are along the bottom edge, tall, and thick, in part to enable comparison to a previous study in Wu2014 (). In order to measure dimensional warping, the prism is removed from the build platform leading to a considerable distortion of the part geometry due to the previously present residual stresses. In Figure 11(a), the experimentally measured displacements in direction, taking on maximal values in the range of , are illustrated. For comparison, in Figures 11(b) and 11(c) the corresponding simulation results either based on a material law with or without plastic hardening, are illustrated. Besides a qualitative agreement of numerical and experimental results, the relative error of the simulation without plastic hardening lies in the range of in the maximal displacement while the variant with plastic hardening overestimates the maximal displacement approximately by a factor of two. For the displacements in direction (not illustrated), the behavior turned out to be (almost) the opposite, i.e. rough agreement of the simulation based on a material law with plastic hardening, large deviations for the variant without plastic hardening. This behavior might indicate that the comparatively simple isotropic constitutive laws applied in Hodge2016 (), and in many other state-of-the-art macroscopic SLM models, are still not accurate enough. Additionally, the simulation results in Figure 12 show the component of the Cauchy stresses resulting from the two chosen preheating temperatures of and . Again, qualitative results can be captured quite well confirming that preheating reduces the peak values of residual stresses.
A slightly different macroscopic SLM model has been proposed in Gusarov2007 (), where the radiation transfer model derived in the authors’ earlier work Gusarov2005 () has been supplemented by a thermal simulation framework based on equations (8). In contrast to Hodge2014 (), no additional phase function field is introduced in order to explicitly solve (8d) for the interface between powder and melt pool. Thus, this interface is only given implicitly by the isothermal line . Because the effective thermal conductivity of loose metallic powders at low temperatures is controlled by gas in the pores Rombouts2005 (), the effective thermal conductivity of the - powders considered in Gusarov2007 () has been chosen to employing the ”gas conductivity” at temperatures below the melting point. Surface contacts between particles are formed when powder melts, so that the thermal conductivity of the powder rapidly approaches the value of the solid material, chosen as with for the investigated stainless steel type 316 L. In the model of Gusarov2007 (), a sharp change from to at the melting point is assumed. Based on the simulated temperature fields and the derived isothermals , the melt pool geometry has been analyzed for different laser beam scan velocities (see Figure 9(a) and 9(c) as well as the corresponding discussion in the paragraph above). In Gusarov2011 (), the pure thermal simulation framework of Gusarov2007 () has been extended to thermo-mechanics incorporating a simple linear-elastic, plain strain version of (10) with elastic constitutive parameters jumping at the phase boundaries.
In Riedlbauer2016 (), the pure thermal problem (8) has been solved for the EBM process based on the finite element method for spatial discretization and an implicit Runge-Kutta method for temporal discretization. Additionally, adaptive mesh refinement strategies have been applied in the direct vicinity of the electron beam. In this model, the penetration of the electron beam into the material is described by means of the model Klassen2014 () and thermal losses due to emission have been considered. The different phases of the investigated Ti-6Al-4V alloy are considered in the thermal model via phase- and temperature-dependence of the thermal material parameters. The heat capacity is prescribed as a nonlinear function of temperature with high gradients around the melting point in order to account for the latent heat of melting (Figure 13(a)). The temperature-dependence of the heat conductivity has been described by a linear law that differs for the three different phases powder, melt and solidified material (Figure 13(b)). In Riedlbauer2014 (), the same authors considered the coupled thermo-mechanical problem (8) and (10), wherein the energy equation (8a) has been derived in a thermodynamically consistent manner from a free energy density functional leading to the well-known thermo-mechanical coupling term of compressible materials not only in the momentum equation (10) but also in the energy equation (8a). The focus of this work lay in comparing the overall computational efficiency of two different numerical coupling schemes, a monolithic coupling scheme solving the thermal and structural mechanics problem simultaneously as well as a staggered partitioned scheme proposed in Armero1992 (). The staggered partitioned scheme allowed for lower computational costs per time step due to the smaller system sizes resulting from a separate treatment of the thermal and structural mechanics problem. However, for stability reasons, the partitioned scheme required a considerable smaller time step size resulting in a higher overall computational effort as compared to the monolithic scheme. Consequently, the authors recommend the latter category of coupling schemes for future application to problem classes such as SLM.
Much work by the group of Stucker et al. Pal2014 () has focused on the development of FEM-based simulation tools for SLM processes, and has been commercialized as the software package 3DSIM. In Patil2015 (); Pal2016 (); Zeng2013 (), the thermal problem (8) is solved, modeling the laser beam as surface heat flux with Gaussian distribution and accounting for heat losses on the powder bed surface due to thermal convection. In this work, an implicit Crank Nicolson scheme for time integration and a spatial finite element discretization combined with a dynamic adaptive mesh refinement and de-refinement algorithm have been employed. While the physical model is comparatively simple and only treats the pure thermal problem, the dynamic mesh adaption algorithm might allow for considerable efficiency gains as compared to uniform meshing strategies. In Zeng2015 (), further steps in improving the aforementioned simulation models are conducted by developing homogenization strategies based on representative volume elements (RVE) for specific part geometries with lattice-like support structures. Among others, these tools have been applied for predicting dimensional warping and suggesting strategies of geometrical compensation.
For the homogenized, macroscopic continuum models as considered in this section, it can generally be expected that high spatial gradients in the thermo-mechanical solution of SLM problems will predominantly occur in the direct vicinity of the laser beam and at spatial locations where also the individual part geometry is characterized by very fine length scales (e.g. transition from high to very low wall thicknesses etc.). Similarly, high temporal gradients are primarily expected in the direct neighborhood of the laser beam position. Consequently, thermo-mechanical simulation approaches based on time- and space-dependent adaption schemes for the spatial (see e.g. Zeng2013 () as discussed above) and for the temporal discretization are very promising to save computational costs.
A further commercial FEM software for SLM simulation is offered by Pan Computing, and was acquired by Autodesk in 2016. For example, the work of Denlinger2014 () extended the contributions discussed in this section so far by allowing for residual stress relaxation. Specifically, the stress and plastic strain of the underlying thermo-elasto-plastic constitutive model were reset to zero, when the temperature exceeded a prescribed annealing temperature. Based on experimental verifications and thermo-mechanical simulations of Ti-6Al-4V parts, it was claimed that temperature-induced stress-relaxation plays a crucial role in EBM processes and that the error in maximal residual stresses and geometrical distortion could reach a value up to as compared to their model without annealing effects. In Denlinger2014a (), the computational efficiency of the model Denlinger2014 () has been increased by employing an inactive element activation strategy as well as a dynamic mesh coarsening algorithm. In the context of FEM-based SLM process simulations, typically two different strategies of material deposition are distinguished: The notion of quiet elements refers to a scheme were all elements representing the final part are already present in the beginning of the SLM process, whereat elements in power layers that have not been deposited yet are characterized by artificial material properties (e.g. very low thermal conductivity , specific heat , Young’s modulus etc.) that are intended to rarely influence the results within the already deposited layers. On the other hand, the notion of inactive elements refers to a scheme that adds new finite elements with every newly deposited layer while elements associated with subsequent layers are not yet present. Of course, quiet element schemes result in a higher computational effort since a large system of equations (containing all degrees of freedom of the final part) has to be solved in every time step and the resulting system matrices might be ill-conditioned as consequence of the artificial material parameters. On the other hand, inactive element schemes are often not supported by commercial FEM codes.
There is an abundant number of further scientific contributions solving the thermal or thermo-mechanical problem prevalent in SLM processes based on macroscopic models and FEM discretization schemes in a manner similar to the examples already given above, with Hodge2014 () being an exception that explicitly tracks the interface between powder and melting phase. While the SLM models discussed so far typically assume temperature- and porosity-dependent thermal properties of the powder based on the initial powder bed porosity, in Cervera1999 () the energy equation (8) as well as an evolution equation for the time- and temperature dependent powder bed porosity field during the sintering process based on a model according to Scherer1976 (); Mackenzie1949 () have been solved simultaneously. These two equations have been coupled by means of a model for porosity-dependent powder bed conductivity proposed in Yagi1957 (). In contrast to this SLS model, an explicit solution of the powder sintering state and porosity in order to accurately determine the current powder conductivity is rather untypical for approaches considering SLM. This may be explained by the fact that the time interval separating the states of loose powder and fully molten liquid is rather short and the modeling error by assuming a simpler ad-hoc relation between powder temperature, porosity and conductivity seems not to be essential compared to the basic assumptions underlying typical macroscopic SLM models. In the works of Dai et al. Dai2001 (); Dai2004 (); Dai2005 (); Dai2006 (), a 3D model for SLS / SLM has been proposed to solve the thermo-mechanical problem (8) and (10) defined by a small number of subsequent powder layers with a special emphasis on multi-material parts. The thermal problem accounts for heat losses due to radiation and convection and employs a Gaussian laser energy distribution. The structural mechanics problem is defined via a constitutive law accounting for elastic, plastic and thermal strains. The effect of volume shrinkage during melting of loose powder is considered by means of a rather heuristic model that moves finite element nodes in direction of gravity by a distance defined via the initial powder porosity and mass conservation. In the work Matsumoto2002 (), a very simple 2D thermo-mechanical model considering one powder layer has been proposed. The mechanical model assumes a linear-elastic state of plane stress and accounts for elastic and thermal strains. The thermal and structural-mechanics problems are solved in a partitioned manner based on the FEM. The contribution Childs2005 () proposes a finite element model for the thermal problem of SLS / SLM processes. There, two different empiric laws have been proposed for the temperature-dependence of powder porosity during the sintering process. Moreover, a porosity-dependent law for the powder conductivity as observed experimentally in Shiomi1999 () has been considered. The laser beam has been modeled as volumetric heat source with an intensity decreasing linearly with the penetration depth. In Roberts2009 (), a computational model for determining the temperature history in a practical SLM process involving the scanning of (up to five) multiple layers has been proposed. The model builds on the previously developed concepts in Cervera1999 () and Kolossov2004 () and extends these approaches to the FEM-based simulation of multiple layers. A comparatively rough yet efficient model is proposed in Zaeh2010 () to estimate residual stresses and dimensional warping effects occurring during the cooling phase of SLM processes. The principal idea is to combine several layers of the SLM process and to apply an equivalent thermal load to this integrated, virtual layer within the heating cycle and to determine residual stresses and dimensional warping effects in a subsequent cooling cycle of the staggered thermo-mechanical scheme. More recent thermal and thermo-mechanical FEM-based SLM models similar to the ones discussed above are e.g. given by Shen2012 (); Shen2012a (); Hussein2013 (); Chen2014 (); Wits2016 (). In Shen2012 (); Shen2012a (), an increased thermal conductivity of the melt pool has been considered to account for the convective melt pool heat transfer, however, without modeling hydrodynamics.
In contrast to the finite element models considered so far, in Verhaeghe2009 () the finite volume method (FVM) has been employed for spatial discretization in combination with an explicit Euler time integration scheme. The underlying model considers the pure thermal problem (8), including the phases of powder, melt, solid and evaporated gas, in order to determine melt pool shapes. Furthermore, the radiation transfer model originally proposed by Gusarov2005 () has been employed. Based on the resulting temperature field, the volume fractions of the phases powder, melt, solid and gas have been determined for each computational cell defined by the FVM. Based on comparisons with experimental data, it was found that the consideration of evaporative heat loss is essential for realistic temperature field results.
An interesting approach that crucially differs from the ones discussed so far has been proposed in jamshidinia2012 (); jamshidinia2013 (). This model somehow combines the procedures typically applied to macroscopic and mesoscopic SLM models. The powder phase is modeled as homogenized continuum with effective, porosity-dependent properties as typical for macroscopic models, whereas the fluid dynamics within the melt pool are explicitly modeled as it is common for mesoscopic SLM models. The melt pool fluid flow is modeled as a compressible, laminar Newtonian liquid flow governed by the Navier-Stokes equations (11) accounting e.g. for temperature-dependent surface tension, temperature-dependent buoyancy forces as well as frictional dissipation in the mushy zone typical for alloys such as the considered Ti-6Al-4V material at the powder-fluid interface. The thermal problem is described on the basis of (8) with a laser beam energy source as proposed in Zaeh2010a () as well as heat losses due to thermal radiation. After solving the thermo-hydrodynamics problem based on a commercial, FVM-based computational fluid dynamics (CFD) solver, the resulting temperature fields are used as input data for a subsequent structural mechanics simulation in the sense of a partitioned field coupling. By employing a commercial FEM code, the structural mechanics problem (10) in combination with a thermo-elasto-plastic material law is solved in order to determine residual stress distributions in the solidified phase. All thermal and mechanical properties relevant for the model have been assumed to be temperature- and porosity-dependent. While the employed laser beam model Zaeh2010a () is only valid for EBM processes, the remaining model constituents seem to be directly transferable to related processes such as SLM. Although the applied field and domain coupling schemes appear to be comparatively rudimentary, the general idea of employing efficient, homogenized macroscale powder models while simultaneously considering the convection-driven heat transfer in the melt pool on the basis of CFD simulations can be regarded as very appealing. Some basic physical phenomena typically resolved by the mesoscopic models discussed in the next section but not accessible by the models considered above, e.g. increased melt pool widths and decreased melt pool peak temperatures as consequence of an outward-directed Marangoni flow, could already be captured by this model, although at a considerably lower resolution - and numerical cost - as compared to the mesoscopic models. A considerable number of similar approaches that might also be relevant for SLM can e.g. be found in the fields of laser and electron beam welding (see e.g. Bachmann2016 (); Chang2015 (); Ganser2016 (); Geiger2009 (); Ki2002 (); Ki2002a (); Li2004 (); Rai2008 (); Semak1999 ()). Future SLM simulation models could strongly benefit from these works.
2.3 Mesoscopic simulation models
While macroscopic models do not consider physical phenomena on the length scale of individual powder particles, mesoscopic models explicitly resolve these scales in order to study melt pool dynamics, melt pool heat transport as well as the wetting of melt on substrate and powder particles. Typically, these models aim at the prediction of part properties such as layer-to-layer adhesion, surface quality and defects on the mesoscale (pores, inclusions etc.). In these models, the initial powder grain distribution is either determined on the basis of contact mechanics simulations, e.g. by employing the discrete element method (DEM, see e.g. Lee2015 ()), or by means of more generic packing algorithms such as rain models Korner2011 (); Meakin1987 (). Subsequently, thermo-mechanical simulations are performed considering the heat transfer within the powder bed, the melting process as well as the heat transfer and the hydrodynamics in the melt pool. Often, the unmelted powder grains are assumed to be spatially fixed and the powder phase is only solved for the pure thermal problem based on a proper laser beam model (see Section 2.1.1 and 2.1.2). The melt pool fluid dynamics are commonly modeled by means of energy equation (8) considering the latent heat of melting as well as the following set of balance equations (11), accounting for conservation of mass and momentum of an incompressible flow:
In the above system of equations, (11a) is the continuity equation requiring conservation of mass, (11b) is the Navier-Stokes momentum equation, (11c) and (11d) represent essential and natural boundary conditions on the boundaries and , while (11f) prescribes the initial velocity field at . The flow is typically assumed to be laminar and incompressible and the fluid to be Newtonian with a stress tensor . Here, represents the rate-of-deformation tensor, is the pressure and is the kinematic viscosity. The body force term typically accounts for gravity. This term might contain additional contributions, e.g. related to the phase change from metal powder to liquid or from temperature-induced buoyancy forces, as discussed below. The equation (11e) accounts for surface tension effects on the free-surface between melt pool and ambient gas. While the velocity field is assumed to be continuous (left side of (11e)) the surface stress exhibits a jump at the interface (right side of (11e)). The jump in interface-normal direction (first term in square brackets) is already present for the case of spatially constant surface tension leading to fluid surfaces with curvature and capillary effects already for hydrostatic problems. The jump in the interface-tangential direction (second term in square brackets) only occurs for spatially varying surface tension values, here solely considered due to spatial temperature gradients, and will for Newtonian fluids (with shear stresses being proportional to velocity gradients) always induce flow. The Marangoni effect resulting from such surface tension variations crucially influences the convective heat transfer within the melt pool. All in all, surface tension significantly determines the melt pool shape and the resulting solidified surfaces (see Khairallah2016 (); Qiu2015 (); Lee2015 ()).
For the numerical solution of the mathematical problem defined by (8) and (11), often a staggered solution scheme relying on a sequential solution of the thermal problem in the powder bed and the fluid dynamics problem in the melt pool is applied. In such a partitioned solution procedure, the powder particle zones with are typically identified as mass increments per time step that are transferred from the solid to the fluid phase due to melting. Consequently, the contour lines define the boundary conditions for the subsequent fluid dynamics solution step. Many of the existing approaches consider a spatial discretization via the finite volume method (FVM) in combination with the volume of fluid (VOF) method in order to track the fluid-gas interface. Temporal discretization of the fluid dynamics problem is predominantly based on explicit time integration schemes. The required resolution of the spatial discretization and the small time step sizes admissible for explicit time stepping schemes lead typically to high computational costs, allowing for an application of these 3D mesoscopic models only to single-/dual-track simulations so far.
Arguably one of the best-known mesoscopic models of this kind has been presented by Khairallah et al. Khairallah2014 (). There, the continuum radiation transfer model Gusarov2005 () as discussed in Section 2.1.1 has been employed as a laser beam source term, while individual particles are resolved on the mesoscale in order to explicitly consider the conductive heat transfer within the powder bed limited by the atmospheric gas in the pores and the point contacts between the powder grains (see Rombouts2005 ()). For spatial discretization, the finite element method based on a uniform Cartesian mesh with element size has been applied. The numerical solution process is based on a staggered solution of the thermal and hydrodynamics problem. While the applied explicit time stepping scheme allowed for a robust treatment of this complex simulation problem in case proper stability requirements are fulfilled, the latter restriction limited the achievable time step sizes and consequently the observable spatial and temporal problem dimensions considerably leading to a computational effort in the range of in order to simulate a single laser track with length dimensions in the range of and time spans in the range of several Khairallah2014 (). The model of Khairallah2014 () includes surface tension effects prevalent in the melt pool flow as well as fluid viscosity and gravity effects, although the latter effect turned out to be negligible given the short time scales and the dominating surface tension effects at the considered length scales. However, the Marangoni effect, i.e., temperature gradient-induced melt flow due to temperature-dependent surface tension characteristics was not considered.
In order to validate the thermal building block of the simulation framework, effective powder conductivities have been calculated with the developed model and compared to analytical solutions and experimental observations. In agreement to Rombouts2005 (), the thermal conductivity of stainless steel 316 powder with a typical packing density of turned out to be only by a factor of three higher than the conductivity of air filling the powder pores. Coupled thermo-hydrodynamics simulations revealed the importance of considering surface tension, which results in a smoother, less granular melt pool geometry and eventually leads to higher effective thermal conductivities within the pool but also from the pool to the underlying substrate. This effect, in turn, led to increased heat transfer, faster cooling rates and consequently smaller melt pool sizes (Figure 14). The rapid thermal expansion directly below the laser beam turned out to induce high melt flow velocities, especially in backward direction. Also the well-known balling effect Levy2003 (); Kruth2007 (); Gusarov2007 (), attributed to the Plateau-Rayleigh instability of a long cylindrical fluid jet breaking up into individual droplets, has been observed for high scanning velocities (see Figure 14(c) and Gusarov2010 () for a theoretical analysis).
In Khairallah2016 (), the model proposed in Khairallah2014 () has been extended to also (implicitly) account for the effects of recoil pressure in the evaporation zone below the laser beam on the subjacent melt pool flow, Marangoni convection, as well as evaporative and radiative surface cooling (see also Khairallah2015 () for more details of the model). Moreover, a more detailed resolution of the laser energy absorption has been achieved based on a ray-tracing model as discussed in Section 2.1.2. The surface tension has been assumed to decrease linearly with temperature, which accounts for Marangoni effects driving the melt flow from the hot laser spot toward the cold rear. As already observed in Gusarov2007 (), the surface temperatures below the laser spot can easily reach boiling values. The resulting vapor recoil pressure adds extra forces to the melt pool surface that create a depression below the laser (Figure 15, bottom right). These observations motivated the differentiation of three different melt pool regions, namely the depression region, which is governed by recoil forces, the tail end region, which is governed by surface tension forces as well as the transition region (see Figure 16).
The approach of Khairallah2016 () does not explicitly resolve vapor flow but employs a traction boundary condition (11d) based on a recoil pressure model proposed in Anisimov1995 (). Accordingly, the recoil pressure depends exponentially on temperature:
In this model, is the ambient pressure, is the evaporation energy per particle, is the Stefan-Boltzmann constant, is the surface temperature of the melting and is the boiling temperature. In order to more accurately model heat losses at the melt pool surface, the evaporation heat flux (again based on Anisimov1995 ()) as well as the heat flux due to thermal radiation emission (see in (9)) has been considered according to
Here, represents a so-called sticking coefficient, is the recoil pressure according to (12), is the molar mass, is the gas constant and again the surface temperature. However, the mass loss due to evaporation has not been considered in the continuity equation (11a) underlying this model. The effect of considering / neglecting these different model components is illustrated in Figure 15. As compared to a continuum radiation transfer model, the employed ray-tracing model can resolve non-uniformities in the radiation absorption across powder particles leading to more narrow, neck-shaped connections to the underlying substrate and in turn to higher heat accumulations within particles as consequence of the reduced thermal conductivity. The ray-tracing model also captures scenarios such as particles that are only partially melted and might contribute to surface and pore defects more accurately. The principal surface tension effect fostering strongly curved melt pool surfaces with minimized surface energy is observable in Figure 15b), where additional melt flow is generated by temperature-induced buoyancy forces. The consideration of temperature-dependent surface tension allows for Marangoni convection, which induces a melt flow from the hot laser spot toward the cold rear (see Figure 15c)). This effect, in turn, increases the melt pool depth, and contributes to melt flow recirculation and melt spattering when colder liquid metal with lower viscosity is ejected from the pool. The additional consideration of recoil pressure due to evaporation leads to a considerably deeper melt pool, increased convective heat transfer within the melt pool and, in combination with the additional heat sinks due to evaporative and radiative surface cooling, to an increased overall cooling rate. Consequently, this case exhibits the lowest amount of stored heat.
In Khairallah2016 (), it has been stated that the overall melt pool dynamics in the direct vicinity of the laser beam is strongly influenced by the recoil pressure, whose magnitude increases exponentially with temperature, leading to the aforementioned depression directly below the laser. The depth of depression due to the recoil pressure is considered to increase with increasing laser power, an effect that is closely related to the so-called keyhole mechanism observed in laser and electron beam welding (see e.g. Rai2008 ()). On the other hand, the fluid dynamics of the cooler back end of the melt pool are governed by surface tension fostering e.g. balling of the tail during cooling Khairallah2014 ().
Via this analysis, Khairallah2016 () enabled also theoretical considerations concerning the creation mechanism of so-called denudation zones. These ”low powder ditches” alongside the melt pool result from lateral particles that are dragged into the melt pool by surface tension and are responsible for the so-called elevated ”edge effect” Craeghs2011 () typically observed for the first laser track of a powder layer. Accordingly, the effects of denudation are stronger in the first track of a layer since powder particles can be attracted from two lateral sides of the melt pool. In Matthews2016 (), this simulation model has been combined with experimental approaches in order to further investigate the physical phenomena and mechanisms responsible for the occurrence of the denudation zone. Accordingly, the ambient gas flow induced by vapor rising above the melt pool may eject surrounding particles from the powder layer and contribute to denudation.
In an alternative model proposed in Lee2015 (), the discrete element method (DEM) is applied in order to study the spreading (i.e., recoating) of Inconel 718 powder accounting for frictional particle-to-particle, particle-to-wall and particle-to-coater contact interaction as well as gravity as driving forces. Subsequently, the thermo-hydrodynamics problem of melt pool (free surface) flow is solved by means of a commercial computational fluid dynamics (CFD) solver based on the finite difference method (FDM), which resolves the fluid-gas phase boundary between melt pool and surrounding atmosphere on the basis of the volume of fluid (VOF) method. Also in this model, temperature-dependent surface tension effects, but no evaporation-related phenomena, have been considered. The laser beam energy absorption within the powder bed is modeled in a simplified manner based on prescribed heat flux boundary conditions on the top surface of the powder layer following the Gaussian distribution of the laser beam. Similar to Khairallah2016 (), a depression of the melt pool directly underneath the laser beam could be observed. Since no evaporation-induced recoil pressure has been considered in Lee2015 (), this observation has essentially been attributed to backward flow of molten metal due to Marangoni convection (comparable to Figure 15c) of Khairallah2016 ()). Also the balling effect has been observed in these numerical simulations and attributed to Plateau-Rayleigh instabilities occurring in dependence of the melt pool dimensions and the powder particle arrangement. Accordingly, a higher packing density was found to reduce the likeliness of the balling effect and to increase the overall surface smoothness while a faster travel speed and lower laser power have been identified as possible sources fostering this effect. On the one hand, Lee2015 () claimed to obtain these simulation results on the basis of a comparable spatial resolution at considerably lower computational effort as compared to alternative approaches such as Khairallah2014 (); Korner2013 (). On the other hand, this reference does not discuss the employed time integrator and the chosen time step size, which have crucial influence on computational performance.
In the review article Megahed2016 (), Megahed et al. discuss several mesoscopic simulation results taken from previous contributions of the authors. The model of Megahed et al. is comparable to the mesoscopic models discussed in this section so far. Among others, it has been shown that lower energy densities yielded smoother melt pool surfaces within the considered range of processing parameters. This effect has been observed in a similar manner in reference Khairallah2016 (), where it has been argued that the melt pool associated with the high energy density reaches higher peak temperature values and, in turn, the surface tension can act for longer time spans leading to a more pronounced balling effect. On the contrary, Lee2015 () has reported about an increasing balling effect with decreasing energy density. The authors of Lee2015 () have argued that the total absorbed energy density decreases when decreasing the laser power at constant velocity, similar, as the total absorbed energy density decreases when increasing the laser velocity at constant laser power, with the latter scenario being known to foster the balling effect and resulting surface roughness. However, the difference between these two cases is that the latter scenario results in a longer, more slender melt pool shape (since less time is available for the melt pool to cool down while the laser beam moves by certain distance along the track), which fosters hydrodynamic instabilities according to the Plateau-Rayleigh theory. All in all, the observations made in Lee2015 () and the observations made in Megahed2016 () and Khairallah2016 () might not necessarily be contradictory. Instead, these observations could be associated with an upper and a lower bound of an energy density interval that allows for stable processing Yadroitsev2007 (); Gusarov2010 (), which is consistent with experimental findings.
In contrast to the continuous laser scan mode as considered so far, some SLM machines use point-wise modulation melting strategies. Megahed et al. have also investigated the melt pool characteristics resulting from such a modulated laser (see Figure 17 for the resulting melt pool shapes at different points in time). It can be seen that the laser modulation and the corresponding exposure time leads to singular deep melt pools. As the melt pools grow in diameter during the exposure time, they join into one continuous melt track on the upper surface of the powder bed down to a depth as typical for a continuous laser beam. However, as argued in Megahed2016 (), there are deeper melt pool peaks that might imply additional anchorage for the new layer to the previous ones. In addition to the discussion in Megahed2016 (), these modulated strategies might also allow for higher energy absorption as the laser beam can penetrate deeply into the pores of the powder layer via multiple reflections instead of directly interacting with the melt pool surface as it is typically the case in continuous scan mode. Potential benefits could arise in terms of lower melt pool surface temperatures and a decreased evaporative mass loss and recoil pressure.
In Qiu2015 (), the melt flow fluid dynamics have been solved by means of the finite volume method (FVM), employing the volume of fluid method (VOF) in order to resolve the interface of the free-surface flow. Besides recoil pressure, temperature-dependent surface tension and temperature-induced buoyancy forces as already considered in Khairallah2016 (), Qiu2015 () additionally supplemented the Navier-Stokes momentum equations (11b) by drag force contributions due to solid/liquid transition in the mushy zone of the melting Ti-6Al-4V alloy based on Darcyâs term for porous media. In the energy equation (8), the additional heat losses due to evaporation, convection as well as radiation emission have been considered. As laser beam heat source, the model proposed in Xu2011 () in the context of laser beam welding has been applied. Since the radiation transfer mechanisms prevalent in SLM might strongly depend on the mode of operation, the application of this, but also of other ad-hoc continuum radiation transfer laser models (see Section 2.1.1), has to be questioned critically. Also in Qiu2015 (), a strong influence of the recoil pressure on the melt pool dynamic in the direct vicinity of the laser beam as well as thereby induced material spattering could be observed, even though significantly decreased laser track sizes as well as structured instead of random powder packings have been considered in this reference. Similar mesoscopic modeling approaches are for example given in Gurtler2013 (); Yu2016 (); Yuan2015 ().
While the mesoscopic models considered so far have been discretized by means of either the FEM, FDM or the FVM, an entirely different approach has been proposed in Korner2011 (). Based on the contributions Attar2009 () and Attar2011 (), the Lattice-Boltzmann method in combination with an explicit time integration scheme has been applied in order to discretize a 2D mesoscopic free-surface, incompressible melt flow model incorporating gravity, surface tension and wetting effects as well as the melting and solidification process prevalent in SLM. The initial powder bed generation is based on a random packing algorithm employing a rain model scheme Meakin1987 (). While the general simulation framework might be applicable to SLM as well as EBM processes, the actual simulations and experiments conducted in this reference focused on the EBM process. Consequently, the energy absorption has been considered to predominantly take place at the powder particle surfaces of first incidence. This fact is reflected by the chosen exponential Lambert-Beer law of light attenuation in dense matter, which is not applicable to the SLM process. The free surface between melt pool and ambient gas is numerically described on the basis of averaged, fluid volume fraction-dependent material properties prevalent in discretization cells at the boundary fluid-gas in a manner similar to the VOM. In Korner2011 (), temperature-dependent surface tension, and thus Marangoni convection, has been neglected. Figure 18 shows the 2D melt pool shape resulting from a spatially fixed laser beam position and different wetting angles , determined by the considered material combination at the triple point solid-fluid-gas, for otherwise identical system parameters.
According to Figure 18, the wetting behavior expressed via the wetting angle considerably influences the resulting melt pool shape. More detailed investigations revealed, however, that the resulting melt pool shape is considerably influenced by the powder packing density and, at least in the range of lower packing densities, also by the individual stochastic realization of different powder beds with identical packing density. Consequently, the local powder topology as well as the wetting behavior resulting from the prevalent material combination solid-liquid-gas might have considerable influence on balling of the solidifying tail. In addition to these single point simulations, also single-track simulations have been performed in Korner2011 (). The melt pool shapes resulting from different scan velocities, once at fixed laser beam power and once at fixed spatial energy density are illustrated in Figure 19. It can be observed that the melt pool stability at constant laser power increases with decreasing laser beam velocity, which, in turn, yields a higher spatial energy density. Small melt pool sizes in the range of the powder grain diameter, which results from low energy input, limits the wetting behavior and fosters round melt pool shapes with minimal surface energy. With increasing size, the relative contribution of surface tension, wetting and gravity effects changes and the melt pool is able to interconnect the powder particles visible in form of a pronounced wetting behavior. On the other hand, comparable melt track topologies are found at constant spatial energy density, advocating this quantity as most relevant factor of influence for the given process conditions. All in all, it has been concluded that the packing density of the powder bed has the most significant effect on the melt pool characteristics.
In Korner2013 (), the model of Korner2011 () has been applied to study a layer-wise buildup process, again for Ti-Al6-V4 powder material. The model simplifies the actual SLM process to 2D by investigating only the plane spanned by the build and scan direction. Furthermore, dimensionless factors characterizing the different physical phenomena and time scales prevalent in SLM have been analyzed in this reference. Based on a comparison of laser beam interaction time and thermal diffusion time, a recommendation for the optimal laser beam velocity has been made in order to avoid overheating and excessive evaporation on the one hand as well as unnecessarily high energy losses on the other hand. It has been argued that surface tension dominates gravity in terms of driving forces whereas inertia effects dominate viscous damping in terms of resistance forces. Furthermore, the Rayleigh time scale, characterizing the relaxation of an interface perturbation under the action of inertia and surface tension forces, turned out to be more than one order of magnitude smaller than thermal diffusion time scales. That is, surface tension-driven coalescence of powder particles follows the melting process nearly instantaneously. Further, in Bauereiss2014 () this model has been employed in order to illustrate the formation of pores. In Markl2013 (); Klassen2014 (), the simulation framework proposed in Korner2011 () has been supplemented by means of alternative models for electron energy dissipation in EBM processes. In reference Ammer2014 () the method has been extended to 3D problems with a special emphasize on code parallelization strategies to increase computational efficiency.
2.4 Microscopic simulation models
The mechanical and physical properties of metals, and therefore of all parts made by SLM, are influenced by their local composition and microstructure characterized by grain size, grain morphology (shape) and grain orientation (texture). In general, the microstructure evolution during solidification processes is governed by the spatial temperature gradients (G), cooling rates and solidification front velocities Kurz1986 (); Glcksman2010 (). As visualized in Figure 20, a qualitative solidification map can approximately be divided into areas of âhighâ and âlowâ cooling rates and into regimes of âcolumnarâ and âequiaxedâ structures. Typically, high cooling rates lead to rapid nucleation of grains and non-equilibrium microstructures, characterized by small grain sizes or small dendrite arm spacings. The competition of thermal gradients and solidification front velocities lead to columnar / elongated, dendritic grain morphologies if high thermal gradients are dominating since in this scenario the liquid-solid interface becomes unstable and preferred orientations grow faster than others to finally form dendritic morphologies. On the contrary, equiaxed, rather spherical grain morphologies are expected if high solidification front velocities are dominating as consequence of a more stable liquid-solid interface.
In contrast to conventional metallurgy such as casting, SLM is based on a highly localized energy supply yielding extremely high heating rates. As consequence of the localized energy supply, the solidified material of previous melt tracks and layers remains at comparatively low temperatures, which in turn yields extreme cooling rates after the laser has passed by. The magnitude of heating and cooling rates typically lies in the range of . The succession of extreme heating and cooling rates eventually also results in high spatial temperature gradients () and solidification front velocities () Majumdar2013 (); Gu2015 (). Such a large range of cooling rates and thermal gradients induces a large variety of microstructures of as-deposited materials differing by grain size, morphology, and orientation. The high cooling rates often lead to meta-stable and non-equilibrium phases and comparatively small grain sizes. As a trend, the grain sizes in SLM are typically smaller and the material strength is higher as compared to alternative manufacturing processes such as casting. With the build direction commonly representing the direction of highest temperature gradients, depending on the employed material often elongated, dendritic grain morphologies accompanied by a higher material strength in build direction can be observed. Additionally, solid state transformations and grain growth contribute to crucial microstructure changes during the repeated thermal cycles (while remaining below the melting point)
experienced by previously deposited material after the initial solidification process (see e.g. the transformation from martensite to alpha and beta phases in Ti-6Al-4V as discussed in Section 3.4).
The knowledge of the resulting micro structure characteristics might provide important details for the formulation of process-specific continuum constitutive laws, which are crucial for accurate residual stress predictions as intended by macroscopic simulation models. On the other hand, global temperature distributions provided by macroscopic continuum models represent essential input variables for studying the solidification process by means of microstructure models. So-called phase field models based on various spatial and temporal discretization strategies can be considered as one of the most established modeling approaches to study solidification processes Wheeler1992 (); Dorr2010 (); Krill2002 (); Chen2002 (); Steinbach2013 (); Apel2014 (); Guo2013 (); Kobayashi1993 (); Qin2010 (); Wang1993 () but also diffusionless solid-solid phase transformations occurring for example during the formation of martensitic non-equilibrium phases Militzer2011 (); Arif2014 ().
Phase field models are based on the definition of a free-energy functional, which typically has the following form:
Here, represents the phase field variable or order parameter taking on the value in the solid phase, in the liquid phase, and a value in the finitely thick interface region of smooth phase transition. Furthermore, represents the (volume / mass) fraction of a certain chemical species prevalent in the problem of interest. For binary systems, the fraction of the second species is given by while for general systems consisting of species, variables are required. As before, is the absolute temperature. Also grain orientation is considered in the phase field model (14). In this context, represents a vector-valued parametrization of grain orientation, e.g. given by rotation vectors, quaternions or Euler parameters. The first term of in (14) represents the temperature-dependent Gibbs free energy either associated with the liquid or solid phase at a certain location in the liquid or solid domain, or a proper interpolation of these values for locations on the phase boundary region. The second term represents energy contributions stemming from liquid-solid phase boundaries characterized by an existing gradient of the phase field. The choice of the parameter determines the resulting phase boundary thickness and requires a compromise between an accurate resolution of the small boundary layer thickness typically prevalent in physical systems (high value of ) and a certain degree of artificially increased thickness for reasons of computational efficiency and robustness (low value of ). Eventually, the last term yields energy contributions due to crystallographic orientation gradients . This contribution fosters uniform growth within individual grains and penalizes the misorientations at grain boundaries, making larger grain sizes with reduced overall grain boundary surface favorable in configurations of thermodynamical equilibrium. Variation of the functional (14) yields the Euler-Lagrange equations of the variational problem, determining the equilibrium solution of the variables :
Starting with a non-equilibrium system, e.g. an undercooled melt, the equilibrium solution defined by the functional (14) is found in practical simulations by transforming (15) into a transient problem and searching for the associated steady state solution. Thereto, the so-called time-evolution phase-field approach applies additional rate terms proportional to , and on the right-hand sides of equations (15). In SLM, the temperature field is required as input variable for the phase field model and might e.g. be provided by a macroscopic model based on (8).
The spatial discretization resolution required for a sufficiently accurate phase boundary representation typically results in a considerable numerical effort. For that reason, many of the available phase field models are applied to 2D scenarios. For illustration of the general capability of phase field models, two existing 3D approaches - even though not applied in the context of SLM - shall briefly be presented. One has been proposed in Krill2002 (), based on a finite difference discretization scheme in space and in time. In this work, only one chemical species (variable not required) and only a discrete set of possible grain orientations is considered. Consequently, the variable representing a continuous spectrum of possible orientations is replaced by additional phase field variables , one for each possible orientation. Figure 21 shows different time steps of a solidification problem allowing for discrete grain orientations. Considering the employed grid and time steps, the computational effort seems to be large. As visualized in Figure 21, the initial number of grains reduces to approximately at the end of the simulation, yielding a decreased overall grain boundary energy in thermodynamical equilibrium.
A further 3D PFM based on a spatial FVM discretization and an implicit backward difference time stepping scheme is given in Dorr2010 (). There, two different chemical species (reflected by one concentration variable ) as well as a continuous grain orientation field as given in (14) has been employed. Figure 22 shows the results of a simulation considering grain growth within a two-phase region (solid phases and ) of a binary alloy (chemical species and ) phase diagram under non-equilibrium conditions. Starting with randomly distributed -phase nuclei (see Figure 22, left), this phase begins growing within the -phase matrix (see Figure 22, middle, for the phase distribution at the end of the simulation). Due to the lower diffusion constant of the -phase and the high cooling rates chosen in this example, the so-called coring effect, i.e. a thermodynamical non-equilibrium state characterized by an inhomogeneous concentration distribution in the -phase (see Figure 22, right), has been observed. The representability of such non-equilibrium configurations is very relevant for the modeling of SLM, given the high cooling rates present. The numerical simulation was based on a finite volume grid and time steps in order to simulate of physical time, leading to computation time on processors.
Only a few approaches to microstructure modeling can be found for SLM/EBM, all of them restricted to 2D. In Gong et al. Gong2015 (), the solidification and growth of primary -phase grains during the EBM processing of Ti-6Al-4V has been studied. Thereto, a phase field model similar to (14) has been employed with representing the concentration of Ti and the concentration of the solute, treated as the combination of Al and V. Thus, the ternary alloy Ti-6Al-4V has been simplified as a binary system, and the influence of grain orientation by means of a parameter has not been considered. The required temperature field has been provided via a finite element solution of the energy equation (8). Figure 23 shows the simulated growth of columnar structures at different time steps during the solidification process. Initially, a pre-defined number of nuclei is placed on the substrate below the molten powder layer. The pre-defined number of nuclei has been determined via a regression equation accounting for an increasing nucleation density with increasing undercooling. Initial dendrite growth occurred primarily in scan direction (horizontal) and build direction (vertical) until neighboring dendrites started to contact each other. Subsequently, vertical growth was observed, resulting in a grain orientation parallel to the build direction as typically observed in experiments.
In Figure 24, the final configurations resulting from different scan speeds are depicted. As expected and in agreement with experiments considered in Gong2015 (), the higher cooling rates induced by higher scan velocities result in a higher nucleation density and, eventually, in a finer grain structure, which is still oriented in the build direction.
A second category of approaches often applied to study solidification processes is given by cellular automaton (CA) schemes. Like PFMs, CA schemes require the solution of the thermal field as input. However, the CA schemes do not rely on an additional phase field variable that explicitly describes the location of the individual phases and phase boundaries. Instead, the evolution of the grain growth is typically constructed in a rather geometrical manner by tracking the solid-liquid interface based on so-called cell capturing schemes under consideration of a given initial nucleus orientation and the current interface geometry. A comprehensive overview of methodologies for the modeling of microstructure evolution, among others considering PFM and CA schemes, is also given by Boettinger2000 ().
In Rai2016 (), the columnar grain growth during solidification of Inconel 718 alloy is studied based on a 2D realization of the CA approach. As visualized in Figure 25, in every time step, the temperature field as well as the phase distribution, are exchanged in a weakly coupled manner. The required temperature field is derived by means of a 2D mesoscopic Lattice Boltzmann simulation model Korner2011 () already discussed in Section 2.3. This study predicts columnar grain growth in the build direction. In contrast to Gong2015 (), different grain orientations are distinguished by this model and the effect of stray grain formation due to partially molten powder particles on the left and right boundary of the computational domain could be taken into account thanks to the powder bed resolution on the mesoscopic scale. Also an epitaxial grain growth in build direction, i.e. grain growth on a newly deposited layer continuing the crystallographic structures of grains prevalent in the previous layer Herzog2016 (), is predicted, leading to grain sizes that considerably exceed the powder layer height. The effect of grain coarsening in the solid phase, i.e. competition between solid grains of different orientation leading to a growth of certain grains, has not been addressed by this CA scheme.
Also in the recent contribution Zhang2013 (), first steps are made towards a 2D simulation framework predicting the resulting microstructure evolution during solidification of 316L in SLM processing. The thermal field has been simulated based on a FEM discretization of (8) while the microstructure model relies on a CA approach.
A problem during the SLM processing of Inconel 718 is the segregation of Nb and the formation of so-called Laves phase grains in the interdendritic region, which may considerably degrade mechanical properties. Nie et al. Nie2014 () employed a numerical scheme based on a stochastic approach for studying the nucleation and growth of dendrites, the segregation of niobium (Nb) and the formation of Laves phase particles during the solidification. Additionally, a FEM-based discretization of the heat conduction equation has been employed for solving the thermal problem, which is required as input for the stochastic microstrucutre evolution model. Based on numerical simulations, Nie et al. demonstrated that low cooling rates as well as a high ratio G/v of temperature gradients (G) to solidification front velocities (v) yield microstructures tending to large columnar dendrite arm spacing and, consequently, continuously distributed coarse Laves phase particles. Thus, it is suggested to increase cooling rates and to decrease temperature gradients in order to obtain small equiaxed dendrite arm spacings and discrete Laves phase particles (Figure 26).
2.5 Summary of existing macro-, meso- and microscale models and potential for future developments
Computational efficiency represents one of the key requirements for SLM simulation tools to be capable of capturing practically relevant time and length scales. Potential improvements of existing simulation frameworks comprise temporally and spatially adaptive time and space discretization schemes, the use of implicit time integrators, efficient code design suitable for high-performance parallel computing as well as the development of efficient and consistent strategies for material deposition and dynamic adaption of the computational domain. Besides computational performance, existing simulation tools require further improvements in terms of model and discretization accuracy. In the following, such improvements but also a potential coupling of the three model classes will be discussed (Figure 27).
Macroscopic models: Macroscopic SLM models typically consider entire parts in order to predict global temperature fields, residual stresses or dimensional warping. The melt pool shape and the resulting temperature distribution crucially depend on an accurate model of laser radiation transfer into and heat conduction transfer within the powder phase. In this context, macroscopic models considering the powder phase in a spatially homogenized sense could be improved by transferring information from mesoscopic models. Also model information related to the melt pool hydrodynamics, which determine the convection-dominated heat transfer within the melting phase, could either be gained from mesoscopic models or explicitly considered in macroscopic models following existing approaches in the fields of laser and electron beam welding. These measures might considerably improve the accuracy of the predicted temperature field in the direct vicinity of the melt pool. While macroscopic models tend to strongly abstract the complex physics prevalent in this regime, they are expected to provide good predictions of the temperature field in regions further away from the melt pool, where heat flow is predominantly governed by the global part geometry, thermal boundary conditions as well as solid material characteristics. However, one of the ultimate goals of macroscopic models is the prediction of residual stress distributions. While temperature-induced thermal strains can be identified as the kinematic origin of these residual stresses, the actual magnitude of these stresses as well as the magnitude of maximally admissible stresses (given e.g. by the yield strength) is essentially determined by the solid material properties, and, thus, by the underlying metallurgic microstructure, which is, however, part of the unknown solution variables of the SLM process. On the contrary, even if derived on the basis of exact temperature solutions, the significance of residual stress predictions is strongly limited as long as comparatively simple material laws are employed. Supplying macroscopic constitutive models with further details concerning material inhomogeneity, anisotropy and temperature-history-dependence, e.g. by means of microscale SLM models, could drastically improve the quality of predictive macroscopic SLM models.
Mesoscopic models: Mesoscopic SLM models resolve length scales of individual powder grains and below in order to accurately describe radiation and heat transfer within the powder bed as well as heat and mass transfer within the melt pool governed by surface tension, wetting effects and evaporation-induced recoil pressure. The ultimate goal pursued by these highly resolved models is to gain an essential understanding of underlying physical phenomena responsible for melt track stability, resulting adhesion between subsequent layers, surface quality as well as creation mechanisms of pores and other types of volume and surface defects. This high spatial resolution naturally goes along with a considerable computational effort. So far, the computational challenge of solving the resulting numerical problems involving complex field and domain couplings has been addressed by explicit time integration limiting the currently observable length scales to single tracks of a few length and the observable time scales to several despite massive usage of high performance computing resources. Concerning model accuracy, recent experimental investigations suggest that the evaporation-induced gas flow above the melt pool might considerably affect melt pool and powder bed dynamics Matthews2016 (). While some first contributions have already considered evaporation in an implicit manner via recoil pressure and heat transfer models, an explicit modeling of the ambient gas flow, at least in the direct vicinity of the melt pool, could allow for considerable insights into the governing physical mechanisms and suggest strategies of influencing these (typically undesirable) gas dynamics. Also effects such as powder particle ejection or denudation are not sufficiently understood and require further investigation. While existing mesoscopic models typically consider spatially fixed powder particles, models that account for powder particle dynamics and collision could allow to study and understand these effects in detail. Eventually, mesoscopic models could benefit from improved thermal and mechanical boundary conditions at the interface between the representative mesoscopic volume and the global SLM part. In this context, macroscopic SLM models might provide useful data.
Microscopic models: Microscopic SLM models investigate the metallurgical microstructure evolution resulting from the high temperature gradients and extreme heating and cooling rates during the SLM process and aim at the prediction of resulting grain sizes, morphologies and textures as well as phase distributions. While mesoscopic models intend rather universal statements concerning optimal adjustment of parameters such as laser beam velocity and power, powder layer thickness and packing density as well as grain shape and size distribution for a given powder material, global temperature distributions, residual stress fields and microstructure evolutions as derived by macroscopic and microscopic SLM models strongly depend on the specific part geometry. Consequently, efficient numerical tools are required in order to enable the simulation of real-size parts in acceptable simulation time and a possible exchange of information between the macroscopic and microscopic scale. Macroscopic models typically provide the temperature field solution as input for microscopic SLM models. Compared to the variety of existing macroscopic and mesoscopic SLM models, microscopic modeling approaches that take into account the specific thermal conditions prevalent in SLM processes can be regarded as being less elaborated. The few existing approaches rely on a simplified 2D representation of the problem although many of the physical mechanisms governing metallurgical microstructure evolution are three-dimensional in nature. Current approaches do not consider more than two individual chemical species. However, an explicit modeling of all relevant alloying elements might be desirable in order to describe experimentally observed effects such as the segregation of individual alloying components. Also mechanisms of solid state phase transformation and grain growth, which have been observed in experiments and are likely to considerably change the material properties during repeated thermal cycles after the initial solidification process, might be desirable supplementations for future models.
3 Experimental studies of thermophysical mechanisms in SLM
In the following, the modeling approaches discussed in Section 2 shall be supplemented by representative methods of experimental characterization. The next Section 3.1 focuses on powder bed radiation and heat transfer as discussed from a rather theoretical point of view in Section 2.1. In the subsequent Sections 3.2, 3.3 and 3.4, the experimental characterization of melt track stability, surface quality and defects, of residual stresses and dimensional warping effects as well as of metallurgical microstructures and grain morphology will be considered. These sections represent the counterpart to the Sections 2.2, 2.3 and 2.4, where exactly these aspects, prevalent at different length scales of the SLM-manufactured part, have been described by means of macroscopic, mesoscopic and microscopic models.
3.1 Optical and thermal characterization of powders
Initial research on powder-laser interaction began not with AM in mind, but rather with surface modification techniques such as cladding and hard-facing processes Haag1996 (). Haag et al. Haag1996 () describe the archetypal experiment to measure absorption coefficients of metal powder, wherein a laser of known power is used to irradiate a powder bed, with the absorption coefficient calculated from the time-derivative of the powder bed temperature. In this study, a laser () was used to heat aluminum, iron, titanium aluminide, and copper powders. Selected particle sizes range from to , validated by SEM. The key results are summarized in Figure 28: In the range of low laser power, material-dependent absorption coefficients as well as the consideration of different powder sizes cause only slight variations in the overall powder bed absorption leading to values of absorption as compared to the incident energy radiation. The increase in absorption of iron and copper at higher laser powers, and correspondingly high temperatures, is attributed to surface oxidation of the material. A practically identical experiment was performed by Tolochko Tolochko2000 () to catalog nominal absorption coefficients for materials relevant to SLS and SLM processing at wavelengths of and . Results indicate that the wavelength represents a superior choice for processing metals, with absorption coefficients generally lying between and , as opposed to values around to for the longer wavelength. Conversely, absorption of radiation is greater for the cataloged polymer materials, and ceramic materials are shown to vary widely.
McVey presented an approach to measure the optical power delivered, reflected, and transmitted through a layer of powder using integrating spheres above and below the powder bed McVey2007 (). In Figure 29, the effective optical depth (the logarithm of the ratio of measured / non-absorbed to incident light), is plotted versus the powder bed thickness for an applied laser beam of wavelength . Accordingly, the absorption of light has a strong dependence on particle size, with finer particles yielding stronger absorption.
More recently, Rombouts et al. Rombouts2005 () have studied the thermal conductivity of powdered 316 stainless steel, iron, and copper. Their apparatus relies on a modulated laser beam to heat an optically opaque sample container, which subsequently conducts heat into a powder bed within the container. Under these circumstances, heat flows through the powder bed to pyro-electric detectors, which are read out using lock-in amplifiers indexed to the laser’s modulation. Experimental results are shown in Figure 30. Accordingly, the powder conductivities are two to three orders of magnitude lower than in the fused material or build platform with bulk conductivities of and for stainless steel, iron, and copper, respectively. It is argued that the thermal transport in powder beds of such fine particles occurs primarily through the gaseous phase (Knudsen diffusion) while direct particle-to-particle conduction is the governing effect for larger particle sizes (millimeter scale).
3.2 Measurement of residual stresses and dimensional warping
As already stated earlier, basically, two thermal regimes responsible for the creation of residual stresses can be distinguished in SLM. First, stresses might be induced in the solid substrate just underneath the melt pool due to high thermal gradients and cooling rates in this region. These residual stresses are influenced by the powder structure and melt pool thermo-hydrodynamics and might result in a delamination of the current material layer. Secondly, residual stresses might also result as consequence of the repeated thermal cycles prevalent in material layers in larger distance from the top layer. The main factors of influence in this region are given by the specific geometry of the part (e.g. slender columns, thin walls and overhanging structures) as well as thermal boundary conditions e.g. in form of fixations on the build platform. Even for primitive geometries such as cubes this second category of thermal stresses might occur due to the temperature gradient in build direction, which typically results in compressive stresses in lower layers induced when the layers of higher temperatures above cool down. In all cases, these thermal stresses need to be reduced, or better avoided, in order to fully exploit the mechanical strength of a part when exposed to external loads, in order to avoid dimensional warping after removing the parts from the build platform and eventually in order to avoid crack propagation typically initiated at locations with residual stress peaks. On the one hand, residual stresses can be relieved by annealing after the SLM process (before build plate removal). On the other hand, such a heat treatment can compromise dimensional accuracy and impose grain growth demanding a tradeoff between these diverse requirements.
Lu et al. Lu2015 () studied the effect of laser scan pattern on the properties of Inconel 718. The laser beam scanning strategy is depicted in Figure 31, where part infill is fused by sequentially scanning square regions known as islands. Specifically, dimensions of , , , and have been studied, using a laser, scan speed, layer thickness, and a spot size. The study begins with a discussion of part density, in which the and specimens are shown to have of their theoretical density. Both the and specimens featured lower densities at and , respectively (see also the next section for a discussion of mesoscale defects such as pores). These findings are congruent with surface micrographs, which show increased porosity within the smaller island sizes, and even cracking within the specimen. Stress-strain data is also included, which shows comparatively little change in offset yield strength and ultimate tensile strength, yet elongation at break increases from to with increasing island size. Again, this is in agreement with the surface micrographs, as the increased density of defects serve as sites to propagate cracks and initiate fracture. Finally, micro-hardness and residual stresses were evaluated. Counterintuitively, the specimen featured the lowest residual stress. This observation is believed to result from stress relaxation due to cracking. As expected, the specimens showed greater residual stresses than the and specimens. The greater thermal gradients of the smaller specimen, arising from a greater degree of heat remaining from the previous pass, can explain the greater residual stresses as compared to the specimens.
In Hodge et al. Hodge2016 (), a combined simulation and experimental approach has been applied in order to access residual stresses as well as the amount of dimensional warping observed in AM components, and to find how those stresses are dependent upon the build orientation and infill scan strategy. As described earlier, Digital Image Correlation (DIC) has been used to measure stress-induced dimensional changes in triangular test components. This experimental technique is based upon time-correlating successive images to extract strain measurements as an object undergoes a physical change, such as in temperature or external loading Khoo2016 (). Hodge uses this technique to monitor the surface of triangular test components as they are removed from the build platform. This process changes the residual stress distribution within the component, thereby inducing a corresponding change in strain. These measurements and previously discussed simulation results are augmented using neutron diffraction experiments to directly measure the residual stress distribution in the test pieces. Like DIC, this experimental method is well established in solid mechanics for mapping stress distributions in crystal-like materials Withers2007 (). The full strain tensor may be determined by approaching the same volume from a plurality of directions. Of course, this measurement represents the average strain tensor over the volume probed, commonly between and . Hodge et al. mainly employed the described experimental techniques for verification of their simulation model (see Section 2.2 for a discussion of the results).
Mercelis and Kruth Mercelis2006 () evaluated residual stresses as dependent on scan strategy, build height, and degree of base plate pre-heating. They relied on sectioning their parts followed by X-Ray diffraction studies to determine residual stresses. X-Ray diffraction determines residual stresses on identical principles to neutron diffraction. However, the comparatively low penetration of the X-Ray technique makes it only useful as a surface technique, hence the need for sectioning. The authors clearly showed that residual stresses increase with component height. As long as the parts stayed connected to the baseplate, very high stress levels, typically in the range of the materialâs yield strength, could be observed. Also the exposure strategy that is being used to fuse the powder layers has a large influence on the residual stress levels being developed. It was found that stresses are larger perpendicular to the scan direction than in scan direction. A scanning strategy subdividing the surface in small islands resulted in a lower maximum stress value. Preheating the baseplate to was shown to reduce induced residual stresses by nearly 10%. Finally, the authors showed that removing the part from the baseplate relaxes approximately 50% of the residual stresses in the part, which negatively impacts the dimensional accuracy of the component once removed.
Havermann et al. Havermann2015 () presented a unique method in which fiber optics were embedded within SLM parts. Thereto, the authors fabricated rectangular test pieces. At the location where a strain measurement was desired, a fiber optic strain sensor (Bragg grating) was placed into the part, thereby enabling the strains induced by residual stresses in the component to also deform the optical fiber. This changes the optical transmission of the embedded Bragg grating, enabling the quantification of strain. The employed small fiber diameter enabled measurements within the first several layers of SLM parts, where interaction with the baseplate is central to the development of stresses. The authors observed increasing compressive residual stresses as the first 6 layers above the fiber were fused. After building 14 layers above the fiber and allowing the component to cool, residual stresses in the range of were measured.
Kempen et al. Kempen2014 () studied the effect of build platform preheating temperature in order to lower thermal gradients and to finally yield crack-free parts with high density. Figure 32 shows three samples, which were built using three different preheating temperatures ( (left), (middle), and (right)). It can be seen that the higher preheating temperature clearly results in lower residual stresses and a lower degree of crack formation. In practical applications, too high preheating temperatures and too long preheating times might however lead to an undesirable sintering of loose powder, which complicates its removal from the final part.
3.3 Characterization of melt track quality and defect detection
In this section, some exemplary experimental approaches for characterizing melt track stability and quality as well as possible defects will be discussed. In recent years, much effort has been invested in monitoring SLM and EBM processes for process and quality control Everton2016 (). Defect detection centers on identifying pores, balling, unfused powder, and cracking as a result of sub-optimal process parameters. In parallel, high speed monitoring has been researched with the aim of real time process control. The results of these efforts have manifested in several commercially available systems for monitoring metal AM processes as shown in Figure 33. In the following, we will discuss several of the key works behind these products, and refer to the review by Everton et al. Everton2016 () for complete treatment.
Craeghs et al. Craeghs2011 () studied melt pool heat transfer using a radiometer and high speed camera, arranged coaxially with the laser. In this work, they describe the so-called edge effect that occurs when a track is scanned without an adjacent contour as already discussed in Section 2.3. Heat from scan segments adjacent to neighboring contours are shown to conduct both into the underlying layers and into the previously solidified, adjacent track in the same layer. However, heat from the initial contour does not have this second path, meaning the bulk of the heat flux must be directed into the underlying layer. As a result, the melt pool is larger, reaches higher peak temperatures, stays hot longer, and incorporates more material than when subsequent contours are melted. Overhanging features were also studied, where it was found that scan parameters suitable for standard processing conditions delivered too much power to the initial layers of overhanging features fabricated above areas of unfused powder. The extra heat enhances Rayleigh instabilities, leading to balling, poor surface finish, and undesirable mechanical properties. Finally, acute corners were shown to present many of the same challenges with temperature control. Specifically, the case of a scan pattern incorporating a 180 degree U-turn was investigated. Again, due to the reduced area for heat to flow out of the melt pool and the partly scan path re-heating, this extra heat results in blob-like formations at the apex of the U-turn. This same experimental hardware was also used to investigate the influence of support strategies on process temperatures for fabrication of overhanging features Krauss2012 (). Here, the authors printed T-shaped structures, in which the overhung regions were supported with uniformly spaced and sized column like additions. The spatial density of the columns was shown to greatly influence the maximum temperatures achieved in the overhanging layer, and therefore porosity and surface finish, as the support structure serves as a heat sink.
Krauss et al. Krauss2014 () performed long-wave infrared (LWIR) imaging to determine component porosity in Inconel-718 specimens. The study begins by deriving an approximation for the thermal diffusivity , which is given by
In this expression, and represent laser power and scan velocity, is the hatch spacing and represents the efficiency of energy absorption in the powder layer (assumed ). Moreover, and are the specific heat capacity and density of the material, as usual. Finally, is the ambient temperature of the build environment, and is the radiometrically determined temperature of the pixel as a function of time . The authors assign an imperfection level based upon the observed porosity and degree of binding error. The imperfection level was determined to be inversely proportional to the diffusivity of cubic specimens at a layer height of with a correlation coefficient of . Based on a proper calibration of the - relation and (16), the proposed procedure enables the in-situ determination of part imperfection based on measured temperature profiles. The method is shown to be sensitive to both delamination and porosity arising from sub-optimal process parameters, with sensitivity demonstrated over a range from nearly fully dense material to a void fraction of 40%. Ultimate sensitivity of this method is not calculated. However, the authors note that the spatial and temporal resolution of imaging sensors represents the greatest challenge. In the follow-up study Krauss2015 (), more detailed information on laser power and scan rate is provided, and the resulting effects on thermal diffusivity, maximal temperature, and observation of melt sputter have been presented.
Simultaneously, Clijsters et al. Clijsters2014 () have studied in-situ monitoring for quality control and are working towards feedback for process control. Using an InGaAs NIR camera and a pair of photo diodes with differing band pass filters, the authors extract signals representing melt pool radiance (intensity of emitted light), area, length, and width. They report success in detection of overheating of the melt pool, specifically at the extremes of fill scan vectors, along with increased porosity in these regions. Using data-fusion techniques, the authors are able to predict the locations of porosity within AM parts, as verified by CT cross sectional images (see Figure 34 for illustration).
Experimentally, process maps can be constituted by gathering data on part density and other properties as related to the scan rate, hatch spacing, laser power, etc. While ultimate control of SLM will arguably require determination of local scan parameters according to the local part geometry and thermal boundary conditions, the present approach is to determine a generalized ”process window” which can be applied to approximately achieve full density. The literature in this field is abundant, considering a wide variety of materials and processes. Thomas et al. Thomas2016 (), propose a methodology of normalizing process parameters for SLM and EBM. The methodology builds on the previous work Ion1992 (), where a dimensionless laser beam power and velocity are defined:
Originally, these dimensionless quantities have been identified from analytical solutions derived for the temperature profiles resulting from a moving heat source on a semi-infinite body, a model that has been applied in order to study heat transport in laser beam welding. In (17a), the effective laser power, defined as laser power times absorption coefficient , is related to the temperature difference between initial powder bed temperature and melting temperature , the thermal conductivity and the laser beam radius . Equation (17b) relates the velocity of the laser beam heat source to the thermal diffusivity . According to Ion1992 (), the choice of the parameters and controls the peak temperature and heating rate. This representation matches Figure 3.
Thomas et al. supplemented the parameters and proposed in Ion1992 () by the following dimensionless quantities:
Here, the dimensionless thickness of the powder layer is defined via the physical powder bed depth and the laser beam radius . Also the hatch spacing is normalized by the laser beam radius. Assuming a powder bed porosity of , the term represents the energy density required to increase the temperature of powder with effective density and heat capacity from the initial temperature to the melt temperature . Again, represents the effectively absorbed laser power. The normalization of the laser power with the product yields an energy per unit volume. By assuming a latent heat of melting in the range of for typically employed metals and preheating temperatures, Thomas postulated the energy density required to increase the powder temperature from to and to melt the powder subsequently as the value . Consequently, the last dimensionless parameter represents the ratio of incident energy density to the energy density required for heating and melting the powder. The parameter combinations found in different experimental contributions from the literature have been plotted in a double-logarithmic diagram as illustrated in Figure 35.
|Heat source (type)||Processing parameters||Thermophysical|
|AM platform||Alloy system||Bed Temp.,|
|Thomas et. al.||Electron Beam|
|Juechter et al.||Electron Beam|
|Vranken et al.||Laser|
|Qui et al.||Laser|
|Xu et al.||Laser|
|SLM 250 HL||Ti–6Al–4V|
|Kamath et al.||Laser|
|SLM Realizer II||316L SS||ASM|
|Unpublished Data||Laser||316L SS||–||–||–||–||–||Al-Bermani|
|Cooper et al.||Laser|
|Boswell et al.||Laser|
|Carter et al.||Laser|