Understanding powder bed fusion additive manufacturing phenomena via numerical simulation

A BSTRACT . The increasing interest in additively manufactured metallic parts from industry has issued a formidable challenge to the academic and scientific world that is asked to design new alloys, optimize process parameters and geometry as well as guarantee the reliability of a new generation of load-bearing components. Unfortunately, understanding the interaction between different phenomena associated to metal-additive manufacturing processes is a very difficult task. In this scenario, numerical modelling emerges as a valid technique to face problems related to the influence of process parameters on metallurgical and mechanical properties of additively manufactured components. This contribution is aimed at summarizing the most important outcomes about metal-additive manufacturing process obtained via numerical simulation with particular reference to powder bed fusion techniques. The fundamentals of additive manufacturing numerical simulation will be also explained in detail. Thermal, metallurgical as well as mechanical aspects are covered.


INTRODUCTION
dditive manufacturing (AM) first emerged in 1987 with stereolithography (SL) from 3D Systems, a process that solidifies thin layers of ultraviolet (UV) light-sensitive liquid polymer using a laser. From then on, different AM processes have been developed so that, in 2010, the American Society for Testing and Materials (ASTM) group A "ASTM F42 -Additive Manufacturing", formulated a set of standards that classify the range of Additive Manufacturing processes into seven categories (Standard Terminology for Additive Manufacturing Technologies, 2012) [1]: VAT Photopolymerisation, Material Jetting, Binder Jetting, Material Extrusion, Powder Bed Fusion, Sheet Lamination, Directed Energy Deposition. Additive manufacturing offers several advantages throughout the design workflow including little lead-time, few constraints, little-skill manufacturing. Furthermore, it is a low environmental impact process because it is characterized by less waste and energy saving. In fact, when compared with traditional manufacturing processes, additive manufacturing can significantly reduce energy usage by using less material and eliminating steps in the production process. Among the different AM technologies, Powder Bed Fusion Processes (PBFPs) are largely used with respect to metallic materials, and for this reason they were chosen as the focus of the present review. PBFPs use a high-power density source to melt or sinter a thin metallic powder layer [2,3]. The solid part is built layer by layer, as schematically shown in Fig. 1. Both selective laser melting (SLM) and electron beam melting (EBM) belong to PBFPs. The first one uses a laser source to melt the powder, the second one uses an electron beam. But compared to SLM, electron beam systems produce lower residual stresses, resulting in less distortion of the printed part and less need for anchors and support structures. Moreover, EBM uses less energy and can consolidate powder layers at a faster rate than SLM, even if the surface finish is typically of lower quality. EBM also requires the parts to be produced in a vacuum and the process can only be used with conductive materials. Despite the numerous advantages offered by PBFPs, the complex interaction between the heat source and the powder layer, as well as the complex thermal phenomena that occurs during the printing process may result in different kinds of defects. Gas pores, for example, may arise from powder surface chemistry modification and/or trapped gas in particles that are released during melting and locked in during solidification. But they may also be due to key-hole effect for deep melt pools. Elongated pores are process induced defects and are due to an inefficient melting regime or spatter and fumes ejection [4]. Unfused power also is a process induced defect while balling [4], that refers to solidification of melted material into spheres, is due to lack of wettability with previous layer, driven by surface tension and directly related to melt pool characteristics [5]. Cracking may also happen because of solidification problems (liquation), residual stresses, surface roughness and other macroscopic defects [6]. Warping can occur between two layers or at the boundary between support and part layer (curling) when the build is stopped and re-started [7]. Delamination -separation between two layers -is another serious defect that cannot be repaired by post processing. It is due to inappropriate melting overlap with previous Heat source Solidified underlying solidified powder or incomplete particles melting [5]. Finally, swelling is a defect similar to the humping phenomenon in welding and it occurs due to surface tension effects related to the melt pool geometry [8].
It is well known that sound parts can be obtained by process parameters optimization. A formidable challenge is understanding the complex interactions between the processing parameters and the metallurgical and mechanical properties of additively manufactured components [9,10]. Process parameters are numerous and collect building direction, layer thickness, platform temperature, spot diameter, source power, scanning strategy, hatch distance, point distance (for non-continuous laser source), scan speed, exposure time. Fig. 2 summarizes some geometrical parameters in spot-to-spot fabrication process that need to be accurately calibrated. For example, although good overlapping of two consecutive tracks is necessary to ensure a good relative density of printed parts, exaggerated overlapping will reduce the printing speed and may be accompanied with negative phenomena like balling effect. On the other hand, the material to be processed is characterized by physical, thermal, metallurgical and mechanical properties that influence the choice of process parameters, as well. Finally, the geometry of the part plays an important role in the fabrication process. Every change in the geometry will change the way that the AM machine performs its fabrication routine affecting the properties of the resulting solid [11,12]. One way to study the interaction between final part properties and processing parameters is to carry out different 'trial and error' experimental tests. This strategy is time intensive, requires extensive raw materials and is not always successful in identifying optimal processing parameters for specific applications. For instance, when a new component (or material) is planned to be produced (or shaped) by additive manufacturing, manufactures print small specimens with varying process parameters until a defects-free specimen is obtained. Results of the tests are often summarized by using the volumetric energy density parameter defined by Eq. (1) [14]: where P is the heat source power, vs is the scan speed, h is the hatch distance and d is the layer thickness. In an alternative formulation d can stand for the beam size. For chosen values of d and h, the volumetric energy density allows summarizing several combinations of power-scan speed values or power-exposure time (with v s = s/exposure time) values in a single 'master diagram' as shown in Fig. 3. Even if the master diagram of Fig. 3 can be used as a guideline in process design, limits in its use are quite evident. It is strongly material and process dependent; so that if the material or a process parameter, different from power or scanning speed, is changed, the master diagram must be change also. Further, as observed by Prashanth et al. [15] and Bertoli et al. [16], the applicability of Eq. (1) is still at stake, even though it has been widely used in literature for optimizing the SLM parameters [17,14]. Eq. (1) may not properly represent the effective energy transferred to melt the powder bed, therefore it needs to be improved involving the material properties and the laser-powder interactions. Another drawback of the above described experimental procedure is that, moving from small samples to real components production, a further improvement of process parameters could be necessary because of the geometric effects that are not taken into account in the parametric study. A second experimental strategy for process parameters optimization is the in-situ monitoring and in-line quality control of the process itself [3,[18][19][20]. However, they are still challenging to implement in actual machines for industrial production.
In the last years, the simulation of AM processes emerged as a powerful tool that allows performing virtual experiments without the need of raw materials and machines [21][22][23][24], sample preparation and extensive characterization methods. For this reason, it might be used to overcome the drawbacks coming from experimental techniques but also to support them by reducing the number of tests necessary for process parameters optimization [25,26]. Finally, it provides a quick feedback loop to designers. It sounds good, if only it were so simple. In fact, modeling all phenomena involved in PBFPs is not an easy task. Analytical models, due to their high time efficiency, may be useful for parametric analyses but are strongly limited to the prediction of thermal fields in very simple geometries. Numerical modeling based on FEA (finite element analysis) or FDA (finite difference analysis) are instead promising methods able to fulfill the need of additively manufactured parts producers, provided that the usability and the computational time of such tools be adapted to industrial applications and not limited to research only. Among the benefits of process simulation, time and cost saving in design of new components, alloy development as well as topology optimization are certainly the most interesting [27,28]. Numerical models were developed in literature with different length scale simulation, microstructural and mesoscale modeling aimed at predicting grain growth and defects [29,30] and macro scale modeling aimed at predicting residual stress and distortion of the as-built parts [31,32]. This contribution is aimed at reviewing different PBFPs simulation works found in literature with the main goal to understanding the influence of process parameters on the fabricated material. The complexity of the process requires a multi-scale modelling approach. Since there is not a common definition in literature of the modeling scale, here we define powder-scale models all those models that simulate the power particles behaviour and the microstructure of one track. In general, stress and strains are not calculated in such 'micromodels', rather, they are used to investigate the effect of powder density and process parameters on microstructure and defects induced by PBFPs. Layer-scale models are those models that are able to capture the thermal and stress fields of one or few layers and may be useful to optimize the scanning strategy or calculate the inherent strain to transfer to the full-scale models. In this case the powder layer is treated as a homogeneous porous media. Finally, full scale models are developed to predict distortions and residual stresses of the real component (Fig. 4). In the following, the paper will be divided into three main sections. The first one is devoted to summarize the major outcomes deriving from powder-scale models; the second one is aimed at describing the layer-scale models and results. Finally, the third part collects the most important numerical strategies adopted to face the computation of full-scale models and summarizes their main outcomes.

Powder deposition
nce the power quality is defined, the powder coating characteristics are the first variables affecting final product quality. By using powder-scale models it is possible to predict the distribution of packing density and microstructure as well as the interaction between the coater arm and the base material [32,33]. The discrete element model (DEM) is used in literature to simulate the powder deposition [34]. The model is based on the solution of the Newton's second law (Eqs. 2 and 3) of motion and rigid body dynamics equation: In Eqs. 2 and 3, m i is the mass of the particle i, is the translational acceleration of particle i, g is the gravity constant, and F ij is the force acting on particle i, given by particle j. For the rotational motion, I i is the moment of inertia of particle i, is the angular acceleration of particle i, ∑ i (r ij x F ij ) represents the total rotational force acting on particle i from the rest of the system. The interaction of particles is calculated basing on Hertzian contact force, F hz [35,36]: with R i and R j radii of particle i and j, respectively, δ the overlap distance of two particles, k the elastic constant, γ the viscoelastic damping constant, nij the unit vector along the line connecting the centers of the two particles, v the component of the relative velocity of the two particles; indices n and t refer to normal and tangential contact respectively. Finally, m eff is the effective mass of two particles. One has to define the particle size distribution, the layer thickness and the rigid box where the particles are inserted. The DEM output is the location and radius of each particle that is transferred to the CFD model as initial conditions. Using DEM and input parameters taken from literature [37,38], O Zhang and Zhang [34] found a packing density of 51% against experimental values in between 50% and 61% [39,40] (Fig.  5). The effects of powder packing densities are studied in Ref. [41]. Authors developed a novel method for the molten pool and porosity modeling in selective laser melting. The general outcome is that, with low packing density values, defects like 'necking' areas may arise ( Fig. 6a) that are detrimental for the laser to penetrate the subsequent powder layer because of the increased layer local thickness. In this situation, lack of fusion defects are induced as schematized in Fig 6b. It is worth noting that a denudation phenomenon can also worsen the local packing density as described in Ref. [42]. According to process parameters, such phenomenon dramatically reduces the powder packing density next to the melt track during the process. This is due to a competition between outward metal vapor flux directed away from the laser spot and entrainment of powder particles in a shear flow of gas driven by a metal vapor jet at the melt track. Finally, Mindt et al. [35] observed that the power size distribution also is an important parameter to be optimized in order to reduce defects in PBF fabricated components. In general, it is found that the lower the mean particle diameter the lower the porosity since the coating process leads to a more uniform powder layer.

Powder-source interaction
Irrespective to welding, in PBF processes the heat source travels over a non-flat surface. The source rays are in fact subjected to shadowing, reflection, diffraction and absorption phenomena induced by the powder layer morphology. This is the reason why in powder-scale models the volumetric heat source is not a good choice. For example, in SLM, heat is generated where the laser beam strikes the particle surface, and then it is diffused into the particle. Even the shadowing effect, where some bottom particles surfaces could be prevented from heating by shadows of other top particles, is a unique feature in PBFPs [34]. By using a volumetric heat source, melting occurs at anywhere inside the particle simultaneously. Thus, neither partial particles melting, nor the shadowing effect could be captured. Li and Tan [43] demonstrated how the light scattering through air voids dramatically change the distribution of energy absorption profile. It is experimentally verified that the total energy absorption of the powder depends on its parameters (size particle distribution, packing density, and so on) and is much higher than that absorbed by the bulk material [44,45]. In literature, the Ray Tracing (RT) Model is used to capture the interaction between the still solid powder layer and the heat source [43,46]. Particles are assumed as rigid spheres fulfilling the relation , were D is the average diameter and λ is the wavelength. The Gaussian power distribution of the beam is then divided in a sufficient number of rays. Each ray is assigned with a certain size, a direction, and amount of power. Reflection on every gas-metal interface and propagation in the gas between two consecutive interfaces is explicitly studied for every ray until it is totally absorbed or reflected outside the power bed. Upon the incidence of a ray on the surface of a particle, a portion of its power is absorbed by the particle surface and the remainder power goes to the reflected ray. The law of reflection determines the direction of the reflected ray, and the absorptivity is calculated by Eq. (5), where γ is the incident angle and n and k are the optical constant of the material. Some authors [34,47] simplified the RT model by considering only the first laser-metal interaction without reflection. Using the RT model, the following conclusions were reached in literature [43,46]: 1) a bed structure has a higher absorption and a deeper laser penetration compared to a flat surface (the powder bed (PB) structure can better enhance the absorption for materials with low absorptivity); 2) absorption distribution along the laser shooting direction doesn't follow an exponential decay, rather, it increases until it reaches a "peak" and drops steeply after that; 3) thin powder bed has higher absorption near the substrate. The total absorption of thin PBs is of the same magnitude as those of thick PBs due to the enhancing effect of the substrate; 4) different laser absorption is calculated using smooth and wrinkled particles models: diffuse mode, due to not mirror-like particle surface, gives the best approximation.

Molten pool
Molten pool and parameters influencing its shape have been studied in literature by fluid-dynamic computation [48]. Vaporization/keyholing, spatter, surface tension-driven fluid flow, recoil pressure and porosity formation are some of the complex phenomena that can be studied by computational fluid dynamics (CFD). Results coming from powder deposition and power-source interaction are used as input for the fluid-dynamic analysis [46,41,34]. The laser heat-driven fluid flow problem is defined by a set of equations such as the conservation of mass (Eq. 6), (ρ is the density, t is the time, u is the flow velocity), the volume fraction conservation (Eqs. 7), (where 1 and 2 are the solid/liquid metal and the gas volume fractions, respectively) and the Navier-Stokes equation (Eqn. 8), (p is the static pressure, μ is the viscosity, ρg is the gravitational body force, ρg[1-β(T-T m )] is the buoyancy force term induced by the temperature dependency of density, K(T) = A mesh (1-f l ) 2 /(f l 3 +ε) is a Darcy condition that suppress the motion of un-melted metal, where fl is the liquid fraction, ε = 0.01 [34] is a small number to avoid zero denominators and A mesh = 10 14 is the mesh zone constant [34]). The last term in Eqn. (8) is the additional momentum source term that in the work by Zhang and Zhang [34] is applied at the metal-air interface [47,49]: In Eqs 9-11 σ is the surface tension, κ is the surface curvature, is the surface normal, is the tangential gradient, is the interface delta function, and is the average density at the interface. P v is the recoil pressure that induces a surface depression over the melt pool. P a is the ambient pressure, λ is the evaporation energy, K B is the Boltzmann constant and T b is the boiling temperature. It is observed from Eqn. (10) that P v is temperature dependent. The temperature values are obtained by solving the time-dependent energy equation in fluid [50]: where (uC p T) is convection within the fluid due to fluid flow, k is the thermal conductivity, Cp is the specific heat and Q is the energy loss because of various reasons (such as convection, radiation, evaporation) and laser heat source [34] represented by the following Gaussian distributed power per unit area [51]: where P is the laser power, α is the absorbance, r0 is the laser radius at 1/e 2 , vlaser is the scan speed of the laser, xi, yi is the initial position of the laser focal center. To simplify the RT approach, in [34] the laser beam power is projected to the top surface of the metal supposing all the heat source is absorbed by the top surface at the first interaction (no optical reflection occurs). Despite the complexity of thermo-fluid dynamic formulation of the problem, interesting phenomena can be captured and studied in detail. For example, Zhang and Zhang [34] showed how the surface morphology of the laser melted power bed is formed by a competition between the Marangoni and recoil pressure effect (Fig. 7). The Marangoni effect is due to the variation of the surface tension (γ) with the temperature. The higher the temperature the lower the surface tension (Fig. 7a). To reduce the surface energy, the hottest fluid is moved to the lower temperature regions. On the other hand, the recoil pressure tends to move the fluid downwards as shown in Fig. 7b. Since the magnitude of the recoil pressure is an exponential function of temperature, whereas the Marangoni convection is caused by a linear temperature dependence of the surface tension, a depression zone is formed in the melted pool. A sufficient good agreement was obtained by Zhang and Zhang [34] between the predicted and observed size of the fusion zone. The higher the scan speed, the lower the melted region depth but with a non-linear trend (Fig. 8). As a matter of fact, it is observed that both the Marangoni effect and the recoil pressure are prone to increase the penetration depth. However, although the energy density is linearly increased, the contribution from recoil pressure and Marangoni force are different because of the different temperature-dependence.  By using fluid dynamic numerical simulation, Willy et al. [46], optimized the overlapping of consecutive molten pool tracks for different laser powers (Fig. 9). It is found that with powers lower than 175 W, a gap is seen between two consecutive molten pool tracks. Zheng et al. [41] used the lattice Boltzmann method (LBM) to solve the Navier-Stock (NS) equations. Compared to finite element method (FEM), finite difference method (FDM) or finite volume method (FVM), LBM describes the fluid dynamics by the collision and streaming of the fluid particles with the great advantage to be much higher efficient than continuum-based approaches (FEM, FVM, FDM). It is found that the surface tension is six orders of magnitude larger than gravity so that fused particles are driven by surface tension rather than gravity to coalesce together (Fig. 10). By using numerical simulation Zheng et al. [41] highlighted the different contribute of interfacial forces on the melted pool morphology. Surface tension, Marangoni convection and recoil pressure effects are shown in Fig. 11. The mere presence of surface tension overestimates the dimension of fusion zone (FZ). When the Marangoni effect is also taken into account, both the FZ size and back material elevation reduce due to the increased convection driven by the tangential Marangoni forces (Fig. 7a). If now the recoil pressure is added to the model, the FZ size decreases further because of the heat loss induced by evaporation while the depression just behind the laser beam and the back-material elevation increase. As depicted in Fig. 12, the laser spot and the depression zone are not found coaxial and this numerical finding was also confirmed by experiments (Fig. 13). For a given laser power, porosities may derive from too high (lack of fusion) or too low (trapped gas) laser speed. Zheng et al. [41] described the lack of fusion porosities formation by numerical simulation (Fig. 14). From 600 mm/s to 1200 mm/s the length of the molten pool increased but the resulting surface profile resulted more irregular. At 1800 mm/s a sudden reduction of the fusion zone dimension was observed, and the melt track became discontinuous. Such discontinuities will result in lack of fusion defects because of the local increase of the layer thickness (as schematized in Fig. 5b). Further increase of laser speed promoted the balling defects formation (  Even porosity induced by trapped gas, known also with the name of keyhole porosity, can be captured with the fluid dynamic model [41]. A too low scanning speed promotes the keyhole regime with a narrower and deeper molten pool (Fig. 15). The recoil pressure becomes much higher than the surface tension so that it could easily cause a disturbance on the surface of the molten pool (Fig. 15a) and entrap gas that will originate a pore defect (Fig. 15c). Finally, it is worth mentioning that, despite the model complexity, the computational efficiency obtained by using the LBM is very high compared to the other powder-scale models found in literature (Tab. 1).

Microstructure evolution
A method used to model the microstructure evolution during solidification of the molten pool is the Cellular Automata (CA). This approach was introduced by Von Neumann in the early 1960s to model complex physical systems [52]. CA is a grid of cells that may be one-dimensional (Fig. l6), two-dimensional, three dimensional and so on. Each cell has a 'state' (for instance, binary state 1 or 0) and neighborhoods. How a cell state changes over time is determined by the state of the neighborhoods. In other words, the cell state, at any moment in time, t, is equal to a function of the neighborhoods states at the previous time, t-1. The CA starts with an initial configuration (IC) (initial state of each cell) and changes over time according to pre-defined rules that will be applied to the lattice each iteration. Thus, the initial state of the lattice will change following the dictate of these rules. In the example of For microstructure evolution simulation, results coming from CFD analysis (temperature and volume fraction) are imported in the CA model, which simulation domain is smaller than the previous one, as a much finer mesh is required to represent the grain structure. In Zhang and Zhang [34], the cell state variable has four possible values: I = 0 (liquid not at interface), I = 1 (liquid cell at the solid-liquid interface), I = 2 (solid at the interface) and I = 3 (solid not at interface). The transformation from liquid at the interface to solid at the interface was modelled using the 'modified decentered square' method [53]. Fig. 18 shows the workflow of the CA method used by Zhang and Zhang [34]. The model was written using Matlab® code while the visualization of results was done by Ovito [54]. Figure 18: Workflow of the CA method with input from CFD results [34].
De Baere et al. [55], used the CA method to simulate the microstructure evolution during a uniform heat treatment at the beta transus temperature of a Ti-6Al-4V additively manufactured sample. In this case the possible cell states were three: 0 = untransformed state (α'), 1 = cell partially transformed, 2 = cell that has finished its transformation (β). The used state transformation rule was: • If any of the cells in the neighborhoods of the current cell has a state "1" change state of the current cell to "1" • If all the cells in the neighborhoods have a state "1", change the status of the current cell to "2" [56]. Fig. 19 shows different historic neighborhoods in two-dimensions grid and that one used in [55] in order to obtain rounder, more realistic grains.  Fig. 20(a) shows the microstructure during the solidification of the fusion zone (grey color). It is found that grains grow from the bottom of the melt pool rather that the left and right of the fusion line. This is because the grain growth rate is proportional to the cooling rate [57], Fig. 20(b). If θ is the angle between the laser scan speed (v) direction and the solidification direction, the solidification rate is given by [34]: Thus, R ≈ 0 at the fusion line at the two sides of the melt pool and R = V at centerline (Fig. 20(c)). The higher the scan speed the lower the grain size due to larger cooling rate. Depending on process parameter the grain structure type can be that one due to 'competitive grain growth' where the grains at the fusion line are oriented in a favorable direction for growth but other grains may surpass the first ones as the fusion line changes its orientation. As a result, the grains at the centreline will grow toward the laser direction. The other grain structure type is called 'centerline grain boundary forming'. In this case grains grow straight from the fusion line to the centerline until they touch each other forming the centerline grain boundary. It is found [34] that by increasing the scan speed, the microstructure changes from 'competitive grain growth' to 'centerline grain boundary forming'.

LAYER-SCALE MODELS
BFP, from a mathematical point of view, can be seen as a transient boundary problem in which the thermal input varies with space and time. It is described by a thermal source traveling over a powder bed, like multi-pass welding [59][60][61]. Filler metal is not present, but the change of properties from powder to continuous material looks like similar to the material addition simulated in welding processes. The temperature histories at each point of the model is then used to predict microstructural evolution according to specific constitutive equations [62][63][64] and finally both thermal and microstructural histories are used as load for the calculation of the induced stress-stain fields [65]. However, it is worth noting that because of the huge number of thin layers and the corresponding heat source runs as well as the fine mesh needed in numerical models to capture the high thermal gradients, when compared to welding simulation the computational cost for additive manufacturing model solving is much higher. In order to reduce the computational time, layer-scale numerical models were developed. Specific AM phenomena like shadowing are not modeled but taken into account using 'ad hoc' material properties and power density distribution functions. In fact, the powder layer is treated as a porous material with isotropic thermal, physical and mechanical properties. The heat source travels over a flat layer. Metal powder deposition is taken into account by particular numerical strategies like activation/deactivation of elements. Advantages are evident. The reduced computational time compared to powder-scale models, allows simulating the complete deposition of one or more layers and obtain information regarding the best scanning strategy aimed at reducing residual stresses and distortions. Such models are desirable for screening and trend prediction purposes, particularly in light of the widening AM processing parameter space. Fluid dynamic may be also taken into account as in the work of Mukherjee et al. [66].

Governing equations
The thermal filed induced by the heat source scanning upon the powder layer is obtained by solving the fundamental equation of heat transfer (2): P (15) where ρ is the material density, Cp is the specific heat capacity of the material (weighted according to the proportions of the various phases), k is the material thermal conductivity (weighted according to the volume fractions of the various phases), T is the temperature, is the temperature rate, L ij (T) is the latent heat (at temperature T) of the i → j transformation, p ij is the phase proportion of phase i which is transformed into phase j in the time unit and that is calculated through phase transformations constitutive equations [64,[67][68] and is a 3D gradient vector operator. The heat transfer boundary conditions of the problem are: q  kT (16) where q is the heat flux at the boundary which, in the PBFP as well as welding process, consists of a prescribed value function of the time and space, convective and radiative heat loss. Thermal and metallurgical results are used as input for the calculation of the stress-strain fields assuming a weakly-coupled relation. The governing equations are: In Eqs. 18-19 [D] is the element stiffness matrix while  e is the elastic term,  th takes into account of the thermal and metallurgical component (volume difference between two phases),  cp is the classic plastic (or viscoplastic) component and  tp is a plastic term induced by phase transformations [69][70][71].

Modeling of materials properties
Convection and evaporation are two heat transfer phenomena that are not negligible in PBFPs simulation [72]. In order to take into account convective heat loss, Mukherjee et al. [66] performed fluid flow calculations in a domain restricted to the molten pool and its adjacent regions to save the total computational time; furthermore, the finest portion of the mesh traveling with the laser source was used to improve the computational efficiency (adaptive mesh refinement). In other works (say, reference [73]), mass convection in the melt pool is approximated by the anisotropically enhanced thermal conductivity method. Since thermo-physical properties of the powder bed depends on alloy and shielding gas properties, size of the powder particles and packing efficiency of the powder bed, effective powder properties are defined as follows [74]: where ρe and Cpe are the effective density and specific heat, respectively; η is the packing efficiency of the powder bed, subscripts s and g stay for solid and gas, respectively.
The effective thermal conductivity may be defined as [75]: where N is the coordination number and L is the ratio of a constant and power particle diameter. The constant value depends on the shielding gas type; for argon, it is equal to 5.4 x 10-4 m-1 while k g is given by the following formula [76]: In [77][78][79] the effective thermal powder conductivity (k e ) is defined as: where  is the porosity (say, 0.4) and k r is the thermal conductivity which is changed by radiation: In Eqn. (25) F is the apparent coefficient, usually 1/3, Dr is particle diameter and B is the Stefan-Boltzmann constant.
In order to simulate the transformation from powder to consolidated alloy, the thermal conductivity is changed when the powder reaches its melting temperature (Tm). In [80] the following relations were used: Ning et al. [81] used a simplified definition for the powder density and thermal conductivity: where γ and β are coefficients that can be taken as 1 [82]. The latent heat can be taken into account by the equivalent capacity method [80,83]: where Tl, Ts and Tv are the liquidus, solidus and vaporization temperatures, respectively, and Lg and Lv are the latent heat (J/K) of fusion and vaporization. In their work, Xiang et al. [84] distinguished also the emissivity of the powder layer from that of the dense material as follows: where ε s is the emissivity of the solid material, A H is the area fraction of the surface that is occupied by the radiationemitting holes and εH is the emissivity of the hole. They depend on powder porosity () according to the following relations [84]: Finally, it is worth noting that some authors [85] modeled the part-powder interface conduction heat transfer by means of a convection heat transfer, eliminating the need for powder elements in the FE model.

Modeling of heat source
Powder-scale numerical models are useful to calculate complex phenomena occurring in the molten pool, but they are steel too time expensive, with some exception [41]. For example, the 4 ms simulation performed in [86] took 140 h. If the molten pool dimensions and the temperature histories at nodes of a numerical model need to be captured in a reasonable computational time, some simplifications are mandatory. Instead of employing laser-ray tracing method in randomly distribute particles, the heat source is assumed as volumetric heat source model. As observed by Zhang et al. [80], heat source models can be divided into two groups, the geometrical modified group (GMG), where the power distribution function try to reproduce the fusion zone shape (like the well-known double ellipsoid volumetric heat source by Goldak et al. [59], used in arc welding simulation (Fig. 21)), and the absorptivity profile group (APG) (applied to Laser PBFP) in which the heat source is described by a two-dimensional Gaussian power distribution on the top surface, while the laser beam is absorbed along the depth of the powder-bed based on the absorptivity function. In any case, each source model needs experiments for parameters calibration. Tab. 2 summarizes the mathematical representation of laser-beam heat sources found in literature [80].  It was found that all the eight heat sources described in Tab. 2 predict molten pool depths significantly smaller than the experimental ones (unless parameters of the heat source models are carefully set by numbers of trial and errors). This was q(x, y,z)  P/(S   OPD  OPD) q(x, y, z)  2 5/2 P  3/2 r l 3 exp 2 attributed to the mass convection in fusion zone that was not taken into account in the models. In order to overcome this drawback, Zhang et al. [80] defined enhanced thermal conductivities as follows: k x   x k k y   y k k z   z k (33) where λx, λy, λz are the anisotropically enhanced factors of thermal conductivity k, which are function of laser power (P) and scanning speed (v): In [80] the following function are proposed: where z and y are the depth and width directions, respectively, and va, a1, b1, a2, b2 are parameters to be calibrated by experiments. Finally, the absorptivity (η) also varies with the combination of laser power and scanning speed used. In [80] the following relation was supposed to be true: where a 3 , b 3 are coefficients to be determined by experiments. Finally, Bruna-Rosso et al. [93] implemented in their numerical model the factitious heat source method in order to take into account the phase change and predict lack of fusion defects as a function of process parameters. That approach, introduced by Rolph & Bathe [94], computes a heat source equivalent to the amount of latent heat either absorbed or released during one-time step at each node undergoing phase change, until completing melting or solidification. Bruna-Rosso et al. [93] used in their numerical model the Goldak's heat source formulation which geometrical parameters (Fig. 21), calibrated trough experiments carried out on 316L stainless steels, were considered functions of laser power (P) and scanning speed (v):

Thermal Results
By using the exponentially decaying equation for the laser heat source modeling and the enhanced thermal conductivity properties (Eqn. 34) and absorptivity (Eqn. 36), Zhang et al. [80] were able to predict the melt pool dimensions with an error ranging from 2.9% and 7.3%. Furthermore, by taking advantage of the Plateau-Rayleigh analysis of the capillarity instability of a liquid circular cylinder [95], it was found that the necessary and sufficient condition of stability is πD/L > 1, where D is the diameter of the cylinder and L is the wavelength of the perturbation. If πD is simplified as the perimeter of the melt pool and L as its length, using results from numerical simulation is it possible to predict the melt pool stability conditions [80] bypassing a fluid-dynamic analysis. Both numerical simulation and experiments showed that lack of fusion defects form during the first laser beam tracks (Fig. 22). This is because a stationary condition between the heat transmitted by the laser beam and the heat absorbed by the layer is not yet reached at the beginning of the process. Powder has a low thermal conductivity and thus the fusion zone is low. As the mass of consolidated material increases, the heat conduction becomes more efficient and the fusion zone increases reaching a stationary dimension. This suggests that a greater heat input would be necessary at the beginning of the process, during the transitory phase, in order to avoid defects. The volume shrinkage induced by the transformation from powder to consolidate material was simulated by different authors [96][97]. The layer is divided into two sub-layers according to the powder porosity. As the powder reached the melting temperature, the upper sub-layer is removed by using for example the 'birth and death' numerical technique. In [84] the volume shrinkage is simulated by using the moving-mesh method. When the temperature at the top surface reaches the melting point and dT/dt > 0, the top surface moves downward with a speed vs that depends on the initial powder porosity. Numerical models that take into account this phenomenon were found more realistic compared to previous ones above all in the prediction of the thermal gradient in vertical direction and molten pool depth [84]. Schwalbach et al. [98] showed in their work how the amount of latent heat be negligible compared to the sensible heat. Using the Stefan number, defined by the ratio of sensible to latent heat: where H sl is the mass specific enthalpy of the solid-liquid transformation and C p is the mass specific heat capacity, it was observed that for Ti-6Al-4V alloy, the value of St = 3.81 means that sensible heat is roughly four times the latent heat of transformation.

Computational efficiency
The computational time depends on the accuracy required, simulated phenomena and numerical strategies used (say, adaptive mesh refinement). Tab. 3 summarized the computational times reached in two layer-scale PBFP models published in recent literature.

Model Description
Computational time Discrete source model [98] Domain: cube with edge length 25 mm 60 s Heat transfer and fluid flow model with traveling fine grid mesh [66] 5 layers, 20 mm long, 5 hatches 5 h Table 3: Computational efficiencies of recent layer-scale published models

Residual stress and distortions
Residual stresses are a self-balanced elastic stress filed caused by a non-uniform elastic-plastic deformation induced by either a thermal gradient (say, welding operations or heat treatments) or an external load (say, rolling, stamping of forging process) [99][100][101]. They are commonly promoted during processes operations and, according to their sign, may severely affect the fatigue resistance and result in geometrical distortions. In PBFP, because of the high thermal gradients and cooling rates (103-108 K/s [102]) induced by the laser or electron beam heat source, residual stress frequently approaches the yield stress of the material [103][104] and may result in severe in-process geometric distortion that may lead to failure of the production process itself [105][106][107][108].
In order to understand how residual stress develops during AM process, the actual deposited layer can be considered first heated to a uniform temperature. At the beginning, the high temperature of the layer promotes compressive thermal stress because its expansion is constrained by the cooler substrate. That compressive stress easily reaches the yield stress of the actual layer because the softening effect induced by the temperature on the consolidated material. When now the layer cools down, its thermal contraction will be constrained by the already cooled substrate promoting tensile residual stress on the just deposited layer. By following this mechanism, it is easy to understand that the in-plane residual stresses are higher than the normal (build direction) stresses. This simply explained mechanism [109] is even much more complicated by the fact that each layer is not uniformly heated/melted. The real non-uniform temperature distribution over the actual layer will depend on the 'scan/raster strategy' and will induce an anisotropic residual stress distribution. Furthermore, the residual stress at each layer will be influenced by the subsequent layer deposition and so on. As expected, the residual stress computation by numerical simulation is highly expensive. This is mainly due to the complexity of the process. An additively manufactured part is in fact obtained by melting hundreds (sometimes thousands) of layers, each one by mans of a laser scanning strategy that can change from layer to layer. The main issue when thermo-mechanical numerical model is to be developed derives from the huge extension of both temporal and length scale. The time step needs to be in the order of micro-seconds against a total simulation time of several hours; the element size should be in the order of μm (to simulate the fusion zone induced by the laser) against the size of the part that is several millimeters. This is the reason why many numerical works in literature focused on simulation of few tracks up to one or few layers deposition (layer-scale models). The main advantages in the use of layerscale models is that no dramatic numerical simplification strategies are required to reach the solution with the possibility to capture interesting phenomena such as the influence of solid phase transformations or the raster strategy on the inplane residual stress. Layer deposition is simulated by activating the elements of the actual powder bed to be scanned while materials properties change from powder to consolidate material when the melting point is reached at elements nodes. Tan et al. [110] developed a layer-scale thermo-metallurgical-mechanical model of Ti6Al4V selective laser melting process. The temperature history at each node is first calculated by simulating the laser scanning over the powder bed. The phase volume fractions (fβ, fα') are then calculated according to the temperature value at each time step by using semiempirical constitutive equations: where β and α' are β-phase and martensitic phase, respectively, B s and B f (1253 K) are the start and finish temperatures of the α to β transformation, Ms (923 K) and Mf are the start and finish temperatures of the β to α' transformation, f'β is the initial volume fraction of the β phase during martensitic transformation. Finally, the effect of specific volume variation induced by solid-state phase transformation is calculated and applied using an equivalent thermal expansion coefficient, which detailed expression is given in [110]: In Eqn. 41 α T is the thermal expansion coefficient, α h and α c are equivalent thermal expansion coefficients taking into account the effect of phase transformation during heating and cooling, respectively. Because of the opposite effect of solid-state phase transformation compared to the thermal deformation, it is found that tensile stress is significantly reduced, and compressive stress increased when phase transformations are modeled. Furthermore, in the four-layer printed part, residual stress parallel to the laser scan direction is negative at the top layer and positive in the rest of the part. Opposite results are obtained by using single-phase alloys . Like in welding, residual stress in the direction of beam motion are found 1.5-2.5 times the magnitude of residual stresses in the transverse direction [115][116][117][118][119]. This is due to the restraint conditions the melt pool is subjected to during heating and cooling (Fig.  23). In the direction orthogonal to the heat source movement the material is free to shrinkage and deform, while in parallel direction the material is constrained by already consolidated material. For the same reason, the in-plane residual stresses are much higher than those in the building direction (orthogonal to the layer) even if the latter tends to increase with the number of layers [120]. The anisotropic feature of the residual stress field can be modified by primary scanning strategy (also referred as the raster strategy) on single and multilayers deposition. This topic was highly investigated in literature by numerical simulation [115,[121][122][123]. With reference to single layer deposition it has been shown that a scan vector reduction reduces the residual stress and distortions [124] while no substantial differences in residual stresses magnitude were found by using unidirectional versus alternating scanning [115]. This phenomenon may be explained by a reduction of the constrain effect with the decrease of the track length, L (Fig. 23) and suggests another raster strategy aimed at reducing residuals tress that is the 'island' scanning method [125][126][127]. Compared to full-length scans, that method has been reported to result in significant reduction in residual stress. Finally, it is worth mentioning that both residual stress and distortions may be reduced by rotating the (full-length) raster by a proper angle layer-to-layer [128]. As residual stresses are caused by thermal gradients, it is reasonable to expect that they are affected by heat source power (HSP) and scan speed (SS), as well. Different authors reported a decreasing of residual stress with decreasing of HSP and increasing of SS [129,130,131]. Similarly, it is found that the higher the energy density the higher the residual stress [125]. This suggests that process parameters should be well calibrated in order to guarantee complete powder fusion and at the same time low residual stress. With regarding to geometry, a thicker base plate reduces the residual stress and distortion because of the increased stiffness of the constraint [132], while the effect of the part size on residual stress was found controversial [116,[132][133][134]]. An increased layer thickness was found to reduce residual stress by reducing the thermal and stress gradient [134,135], Powder Solids Restrained Restrained Laser direction Free L while a wide number of studies confirmed by FE modeling, agree in stating that the base plate preheating dramatically overshadows all other process parameters. In particular, it is found that the higher the base plate temperature, the lower the residual stress, with reductions ranging from 10% to 40% [136][137][138][139]. The higher base plate temperature reduces the cooling rate that in turn reduces the residual stress; however, it also increases the grain size and may compromise the mechanical properties according to the Hall-Petch relationship. This is also the reason why parts made out of electron beam melting (EBM) are characterized by lower values of residual stress compared to SLM components. As a matter of fact, EBM allows using higher preheating temperatures in the range of 300-1100 °C compared to 100-400 °C of SLM processing. Another way to reduce the cooling rate is to pre-and/or re-scanning the melted layer. Using this technique, the literature reports a decreasing of residual stress from 6% to 55% [126,140,141]. The interlayer dwell time is another parameter affecting RS with changes of up to 55% when adding a 40 s dwell [142]. This parameter was however found to be heavily material dependent producing opposite results for IN625 and Ti-6Al-4V.

Beam variables
Power [129][130][131] The lower the power the lower the residual stress Volumetric Energy density [125] The lower the volumetric energy density the lower the residual stress Scan speed [130,131] The higher the scan speed the lower the residual stress Scan strategy + The XY alternating strategy induces the lowest and the most uniform residual stress distribution [115,[121][122][123][125][126][127]128]

Process conditions
Among process conditions, pre-heating temperature is the most effective in reducing residual stress [138] Material properties Effect of solid-state phase transformation on residual stress [110]  The higher the diffusivity the lower the residual stress  The higher the thermal conductivity the lower the residual stress  The lower the yield stress the lower the residual stress [143]  The material properties affect the residual stress field indirectly by influencing the thermal filed through the thermal diffusivity and, directly by the variation of the yield stress with the temperature. Despite a high dispersion in results due to the variation of process parameters and geometries, Bartlett and Li [143], by analyzing a lot of literature works, found a strong linear relationship between the residual stress and the material thermal diffusivity and conductivity. In particular, the higher the diffusivity the lower the residual stress (according to the solidification theory). A weaker linear relationship was also reported between RS magnitude and UTS and YS, where increasing material strength generally resulted in increasing residual stress. It is worth mentioning that a high coupling is observed between the material properties and the process parameters, so that universal relations between process parameters and residual stress are not possible. For example, Denlinger et al. [142] studied the effect of interlayer dwell time on residual stress of SLM parts built out of both Inconel 625 and Ti-6Al-4V. They found completely opposing responses to increasing the dwell time, where RS decreased in IN625 and increased in Ti-6Al-4V. Chemical composition of the alloy may affect the residual stress, as well, by the effect of solid-state phase transformations [110,144] or by choosing an ad hoc eutectic composition that keeps the material in a semi-solid state at high temperature [145]. Fig. 24 collects the most important outcomes about the influence of process parameters and material properties on residual stress induced by PBFP.
It has to point out that most of the works in the literature assumes a weakly-coupled (instead of a fully-coupled) relation to calculate residual stress [146,147]. First thermal history at each node is calculated and then applied as load for the mechanical analysis. Layer-scale models [121,124,137,148] are almost mature in their form but their validation is often obtained by measuring distortion of the printed parts and less frequently by experimental residual stress coming from other works that are still rare in literature. Finally, Moser et al. [149] outlined how temperature dependent thermal material properties do not influence significantly the residual stress field [150,151]. Rather, like in welding simulation, the most important parameter affecting the residual stress value is the temperature dependent yield stress value. This highlights the challenges of making credible predictions in PBFP models.

FULL-SCALE MODELS
t is quite easy to guess that solving a full-scale model by using the layer-scale numerical technique is quite impossible with the actual computational power, so far. Therefore, deep simplifications of the thermo-structural problem are developed in literature. Using the dynamic mesh coarsening the number of nodes required to model the whole printed part is reduced since at each time step the finer mesh zone in the surrounding of the laser spot moves with the laser itself [142]. In this way, the accuracy of standard meso-scale models is maintained but the fine time discretization remains. To date, the most efficient strategies used to solve full-scale models are:  inherent strain-based approaches, where a plastic strain filed calculated by experiments or a layer-scale model is applied layer-by-layer on to a macro-scale simulation of the whole part;  agglomeration (or lumping) approaches, where many layers are lumped into one larger computational layer (or block) and the thermo-mechanical computation is carried out in the full-scale model like in layer-scale model. The inherent strain method was first introduced by Ueda et al. [152] in 1975 for rapid solution of welding numerical models [153]. As schematized in Fig. 25, the inherent strain is the plastic permanent strain the causes residuals stresses. In other words, by imposing the inherent strain as plastic strain tensor on the part and by imposing the equilibrium conditions it is possible to obtain the residual stress field with a linear elastic analysis. With application to AM process simulation, this approach involves calculating the inherent strain using experiments [154] or a meso-scale model with a sufficient number of layers in order to reach the steady-state conditions. The inherent strain is then applied as an initial condition over a full layer on a macro-scale simulation of the full part, layer-by-layer according to the building strategy, with the advantage to carry out at each layer activation a simple liner elastic analysis [155], thus bypassing the thermal analysis. The limitation of the inherent strain method is that it doesn't take into account the heating up of the part during the build; different thicknesses of the same printed part experiment different thermal histories and therefore different plastic strain fields [156]. Keller and Ploshikhin [157] found a good agreement between experimental and numerical distortion of the part obtained with the inherent strain method. The inherent strain was calculated by simulating the laser scan of only one layer after having calibrated the laser heat source parameters by experiments. The different raster strategies were obtained by rotating the inherent strain tensor at each macroscopic layer according to the building strategy (Fig. 26). I Figure 25: Inherent strain definition: a) without residual stress, b) with residual stress; c) removal of residual stresses by cutting the part. ε* is the inherent strain. Figure 26: Comparison between experimental and numerical distortion of a cantilever specimen obtained by using the inherent strain method [158].
The second kind of approach used to solve the full-scale model is the thermo-mechanical simulation of the entire part but using a layer agglomeration strategy to speed up the calculation. Each macro layer (or block), when activated is either heated at one time [158,159] or scanned by a heat source with fictious parameters in order to reach the full penetration of the block itself [160,161]. A 'multiscale approach' is also adopted in literature [147,162] where the equivalent body heat flux is first derived from a micro-scale laser-scanning model and input to a meso-scale model for the simulation of the in-plane stress and finally into the full-scale model to obtained the final residual stress state. Li et al. [163] looked at the effect of macro layer thickness on the residual stress and distortion. They found that a scale-up of the layer thickness from 30 μm to 1.5 mm dramatically reduces the computational cost by more than 3 times compared to a 0.75 mm macro layer thickness. They confirmed that orthogonal scanning pattern between two adjacent block layers is beneficial for reducing stress and distortion. In general, thermo-mechanical computation of the whole part, even when employing an agglomeration approach, is more time consuming compared to the inherent strain method but even more reliable since the effects of part geometry on the thermal history are at least partially taken into account. a b c CONCLUSIONS umerical modeling of additive manufacturing processes has reached a good maturity most likely thanks to the past experience acquired in developing welding numerical models. The complex phenomena involved in additive manufacturing process, such as powder deposition, shadowing, material vaporization, powder consolidation, thermal fields, residual stress and distortions require a multi-scale modeling. When focusing on process parameters inducing defect and microstructure, powder scale models are of great help as they are able to predict the conditions causing defects like lack of fusion or keyhole porosity. If the influence of process parameters on residual stress and distortion of the printed part is of primary importance, layer-scale and full-scale models are the best choice. Numerical modeling has the advantage that is possible to virtually study the effect of one process parameter on defects, microstructure and residual stress by keeping constant all the other ones. It was the main objective of the present review that wanted to focus on the main outcomes coming from numerical models of powder bed fusion processes found in literature. To date, the simulation cannot substitute completely process parameters experimental calibration but allows reducing dramatically the number of trials. The main challenge is to enhance the computational efficiency of the numerical models, whatever the scale of investigation.