Physics-Based Modeling and Predictive Simulation of Powder Bed Fusion Additive Manufacturing Across Length Scales

Powder bed fusion additive manufacturing (PBFAM) of metals has the potential to enable new paradigms of product design, manufacturing and supply chains while accelerating the realization of new technologies in the medical, aerospace, and other industries. Currently, wider adoption of PBFAM is held back by difficulty in part qualification, high production costs and low production rates, as extensive process tuning, post-processing, and inspection are required before a final part can be produced and deployed. Physics-based modeling and predictive simulation of PBFAM offers the potential to advance fundamental understanding of physical mechanisms that initiate process instabilities and cause defects. In turn, these insights can help link process and feedstock parameters with resulting part and material properties, thereby predicting optimal processing conditions and inspiring the development of improved processing hardware, strategies and materials. This work presents recent developments of our research team in the modeling of metal PBFAM processes spanning length scales, namely mesoscale powder modeling, mesoscale melt pool modeling, macroscale thermo-solid-mechanical modeling and microstructure modeling. Ongoing work in experimental validation of these models is also summarized. In conclusion, we discuss the interplay of these individual submodels within an integrated overall modeling approach, along with future research directions.

The physics-based modeling and predictive simulation of these processes offers the potential to advance fundamental understanding of the governing process physics on different length scales and mechanisms of process instability and defect creation, which link process and feedstock parameters with resulting part and material properties.
The present article gives an overview of recent research activities and developments in these fields by our research group at the Institute for Computational Mechanics of the Technical University of Munich (TUM) together with our collaborators from the Massachusetts Institute of Technology (MIT), Lawrence Livermore National Laboratory (LLNL), the University of Innsbruck, and the Hamburg University of Technology. Beyond contributions to the individual disciplines of partscale thermomechanical modeling (Section 4), mesoscale melt pool (Section 3) and powder modeling (Section 2), and microstructure modeling (Section 5), the question will be addressed how information between these different models can be exchanged and individual insights can be combined to achieve a holistic picture of the entire process chain in metal PBFAM (Section 6). All models have been implemented in our in-house parallel multi-physics research code BACI [63] jointly developed at the Institute for Computational Mechanics. Also the corresponding simulation results have been produced with BACI unless stated otherwise.

MESOSCALE POWDER MODELING
The characteristics of feedstock powders, and the resulting spreading kinematics critically influence the packing density and quality of the powder layer in PBFAM, which in turn is coupled to the melting process [64,65]. In the authors' recent contribution [66], a model for cohesive metal powders has been proposed and was subsequently applied to study the powder spreading process in metal PBFAM [29,67]. Important model equations and exemplary results are recapitulated in the following.

Model equations
The proposed model is based on the discrete element method (DEM) [68][69][70][71] and represents bulk powder mechanics on the level of individual powder particles in a Lagrangian manner. This model describes the kinematics of powder particles in 3D space by six degrees of freedom, i.e., the position vector of the particle centroid as well as the rotation vector and associated angular velocity . In the present work, the powder particles are assumed to be spherical, which is a good approximation for the considered plasma-atomized powders, and follow a log-normal type size distribution [66]. Of course, a DEM representation of more general particle shapes, e.g., based on multi-sphere approaches, is possible as well. The equations of motion of a particle follow from the balance of linear and angular momentum: Here, = 4∕3 3 is the particle mass, = 0.4 2 is the moment of inertia of mass with respect to the particle centroid , is the particle radius, is the mass density, and is the gravitational acceleration. The interaction forces and torques between particles and on the right-hand side of (1) consist of normal contact forces , tangential (frictional) contact forces , adhesive forces , and rolling resistance torques . Here, ∶= − is the vector from the center of particle to the contact point with particle . As shown in [66] and the subsequent work [29], adhesive forces, mainly resulting from short range van der Waals (vdW) interactions between particles, dominate the dynamics of micron-scale metal powders and critically drive their spreading behavior during the powder recoating process in PBFAM. In [66], an inter-particle force law, following the Derjaguin-Muller-Toporov (DMT) model, has been proposed: Based on (2), the adhesive forces read = ( ) , with the normal vector = ( − )∕|| − ||. Moreover, is the surface energy, is the Hamaker constant, ∶= + is the effective interaction radius, 0 is the distance at which the vdW force equals the pull-off force 0 , and * is defined as the cut-off radius at which the vdW has a relative decline of 0 ∶= ( * )∕ 0 with respect to the pull-off force 0 (taken as 0 = 1% in our studies). In [66], the surface energy has been identified as most critical model parameter because: (i) it might vary by several orders of magnitude for a given material due to variations in particle surface roughness/topology and chemistry, and therefore cannot be taken from standard data bases but rather has to be identified/calibrated for the specific powder material; and (ii) it uniquely defines the magnitude of the particle pull-off force 0 , which is critical for the flowability of these powders. For effective surface energy values in the range of = 0.1 ∕ 2 as identified for these powders [66] and mean particle diameters in the range of = 2 ≈ 30 , the adhesive forces according to (2) are already by one order of magnitude higher than the gravity forces acting on these particles. Figure 2 represents different powder rheological setups along with their experimental and model-based realization. Specifically, funnel tests according to Figure 2 (a) allow for static angle of repose (AOR) measurements, rotating drum tests according to Figure 2 (b) allow for dynamic angle of repose (AOR) measurements, and rotating-vane rheometers allow to measure reaction torque curves as function of angular velocity Ω and compressing normal force , which both can be controlled in these experiments. The purpose of such rheometer studies in the context of metal AM powder spreading process is two-fold: (i) to  identify, i.e., calibrate, unknown parameters underlying the DEM model, which is essential to capture the physical powder behavior with sufficient accuracy; and (ii) to determine if the spreadability of a given powder material can be directly predicted by measurement of key rheological parameters such as static/dynamic AOR or powder characteristics gained from the rotatingvane rheometry. As emphasized in the last section, specifically the precise calibration of the effective surface energy associated with the given powder material is of highest importance to achieve an accurate and predictive powder model. The simulation results presented in the subsequent section are based on a calibrated model, whose effective surface energy has been identified from experimental static AOR measurements resulting in a value of = 0.1 ∕ 2 [66].

Powder Spreading
In [29], a powder model as described above has been employed to study the powder spreading process in metal PBFAM using a rigid blade as recoating tool (see Figure 3 (a)). The mean values <...> as well as the standard deviation (...) of the spatial packing fraction field Φ( , ) and the spatial surface profile ( , ) have been defined as metrics (see [29] for the exact definition) to rigorously assess powder layer quality. Special focus was on the influence of powder cohesiveness, representing powders with different mean particle size, on the resulting layer characteristics. Typical results are shown in Figures 4 (a), 5 (a) and 5 (b). Accordingly, the powder layer quality decreases with increasing cohesiveness, which becomes visible through decreasing mean values of the packing fraction field and increasing standard deviations of packing fraction and surface profile, with the latter being a metric for surface roughness. In addition, the influence of the nominal layer thickness 0 as multiple of the (theoretical) maximal particle diameter ,0 (e.g., taken as ,0 = 50 for a powder size distribution with 90th percentile 90 = 44 ) is illustrated in these figures. For example, the mean packing fraction considerably increases with increasing nominal layer thickness approaching a saturation value at approximately 0 ∕ ,0 ≈ 4. Figure 4 (b) shows recent results, where packing fraction values from spreading simulations and experimental spreading studies are compared, employing a novel X-ray microscopy technique for packing fraction measurement. Accordingly, experiments and simulations are in very good agreement, confirming that packing fraction increases with increasing nominal layer thickness 0 and approaching a saturation value in the range of 0 ∕ ,0 ≈ 4. Note that slightly different definitions of the reference volume underlying the packing fraction calculations in Figures 4 (a) and 4 (b) have been applied. In addition to the rigid blade studies discussed so far, Figures 3 (b) and 3 (c) show results from recent spreading studies employing also flexible blades and rotating rollers. Specifically, the flexible blade simulations were enabled by a staggered coupling algorithm for structure-particle interaction [72] taking advantage of the finite element method (FEM) and discrete element method (DEM) modules in our research code BACI. In Table 1 , first results are depicted for the mean values and standard deviations of packing fraction and surface profile resulting from simulations with rigid and flexible blade respectively non-rotating and counter-rotating roller.
The simulations have been carried out for the most cohesive powder depicted in Figures 4 (a), 5 (a) and 5 (b) (i.e., ∕ 0 = 4) with ,0 = 50 and 0 = 3 ,0 . To predict an upper limit for the maximal achievable layer quality, particle-substrate  (a) Influence of powder cohesiveness and nominal layer thickness on resulting mean value of spatial packing fraction field Φ( , ) [29].
(b) Experimental measurement and associatd computational model prediction for mean value of spatial packing fraction field Φ( , ) [67].

FIGURE 4
Mean value of spatial packing fraction field Φ( , ) for different powders and layer thicknesses.
adhesion has been set identical to particle-particle adhesion and particle-blade/roller adhesion has been set to zero [29]. As expected, the mechanical compression forces induced by the roller geometry lead to a considerably improved layer quality as compared to the rigid blade. Interestingly, the results from the non-rotating and counter-rotating roller are almost identical. This might be explained by the idealized wall-adhesion properties in this study. With non-zero particle-roller adhesion and potentially decreased particle-substrate adhesion, which might yield considerable stick-slip motions of the powder on the substrate [29], a stronger dependence of the layer characteristics on the roller rotation is expected. Also the low layer quality resulting from the flexible blade is an unexpected result at first glance. However, from looking at the dynamic spreading motion it would become apparent that the powder-blade interaction leads to dynamic bending oscillations of the flexible blade, which explains the irregularity of the resulting powder layer. At this point, it is emphasized that the flexible blade setup considered here should rather serve as a proof of principle for the proposed modeling and simulation framework, while flexible blades applied in practice are typically of a geometrically more complex shape.  Figure 4 (b) demonstrates experimental validation of the DEM powder simulations using an X-ray microscopy technique recently developed at MIT [67]. This method was developed to provide direct interrogation of local powder deposition, in contrast with optical interrogation techniques (e.g., [73][74][75][76][77]) that accurately assess layer topography but cannot reliably determine the volume of material deposited as packing density is also spatially varying [77,78]. Herein, etched silicon templates, manufactured to sub-micron tolerances of nominal depth, flatness, and surface roughness, provide precision control of powder geometry. Specimen powder layers are created by sweeping powder into the template, and any variation in powder deposition may be fully ascribed to stochastic powder flow and not to disturbances from poorly controlled boundary conditions. X-ray transmission of the powder layer is assessed by first imaging the empty template and again after creating a layer specimen; these data are interpreted as an effective thickness of metal powder using a radiation transport model considering: the polychromatic X-ray source emission spectrum, wavelength-dependent absorption spectrum of objects (e.g., template and powder layer) in the beam path, and detector physics including scintillation and signal gain. Using this technique, experimental conditions may be matched as closely as possible to simulations for validation. The aforementioned figure demonstrated close agreement in which commercially obtained  Ti-6Al-4V powder has been spread into layers nominally 50 to 250 thick. After careful model calibration via static angle of repose as specified above, we observe that both the experimental measurement and model prediction of packing fraction asymptotically increase with increasing layer thickness as particles are able to assume more dense configurations as compared to cases where layer thickness approaches (particle size). More broadly, the experiment reproduces key relationships between powder flowability, where high angle of repose powders create comparatively sparse layers with high variance in deposition.

MESOSCALE MODELING OF MELT POOL THERMO-HYDRODYNAMICS
In the following, the most important model equations and exemplary results of a novel smoothed particle hydrodynamics (SPH) formulation for thermo-capillary phase change problems will be recapitulated, which has recently been proposed by the authors and applied to metal AM melt pool modeling [25].

Model equations
The domain of melt pool thermo-hydrodynamics in metal AM is split into the liquid melt phase Ω , the atmospheric gas phase Ω as well as the solid phase Ω , allowing for reversible phase transition between liquid and solid phase. The liquid and gas phase are governed by the weakly compressible, anisothermal Navier-Stokes equations in the fluid domain Ω = Ω ∪ Ω . The problem is described by a set of six equations, namely the continuity equation (3a), the Navier-Stokes momentum equation (3b) (3 components), the energy equation (3c) and an equation of state (3d) associated with the weakly compressible approach.
Note, that the energy equation is solved for the entire domain Ω = Ω ∪ Ω ∪ Ω , i.e., also for the solid phase. The six primary unknowns are given by velocity (three components), density , pressure and temperature . The momentum equation (3b) contains contributions from viscous forces = ∇ 2 with dynamic viscosity , surface tension forces̃ , wetting forces̃ , evaporation-induced recoil pressure forces̃ , and body forces . Here, the superscripts of interface forces refer to the corresponding interface, e.g., to the liquid-gas interface (lg) or the triple line solid-liquid-gas (slg). In particular, evaporation is modeled on basis of a phenomenological model for the evaporation-induced recoil pressure as proposed by Anisimov [79]: with the atmospheric pressure , the molar latent heat of evaporationh , the molar gas constant , and the boiling temperature . Here, is the normal vector of the liquid-gas interface and the corresponding surface delta function distributing interface surface forces across a diffuse interface domain of finite thickness. The energy equation (3c) contains the mass-specific heat capacity , the heat flux ∶= − according to Fourier's law (with thermal conductivity ) as well as heat fluxes stemming from the laser beam heat sourcẽ and from evaporation-induced heat losses̃ . The former is modeled as a surface heat flux with Gaussian profile. The latter is consistent with the aforementioned recoil pressure model [79] and reads: where the enthalpy rate per unit area results from the vapor mass flow per unit areȧ at the melt pool surface and the sum of the specific enthalpy ℎ( ) and the latent heat of evaporation ℎ , both per unit mass. Moreover, ℎ,0 is a reference temperature of the specific enthalpy and is the molar mass. Finally, ( ) is the recoil pressure defined in (4) and the so-called sticking constant which takes on a value close to one, i.e., ≈ 1, for metals [19,22]. Finally, in the equation of state (3d) the reference density 0 , the reference pressure 0 = 0 2 and the artificial speed of sound can be identified. Note that the commonly applied weakly compressible approach only represents deviations from the reference pressure, i.e., 0 = 0, and not the total pressure. Equations (3) are discretized in space by smoothed particle hydrodynamics (SPH) and in time by an explicit velocity-Verlet scheme. It is emphasized that SPH is a Lagrangian approach, i.e., material particles are used for discretization that are convected by the fluid velocity and directly carry the phase information without requiring additional phase tracking schemes. In particular, after spatial discretization the model equations (3) (in strong form) are evaluated at the positions of these discretization particles taking into account the material parameters corresponding to the respective phase of a particle, i.e., solid, liquid or gas phase. Continuous primary variable fields are approximated by applying a smoothing kernel function to the discrete particle values of these variables, which is required to define spatial gradients.

Exemplary simulation results
While the accuracy of the individual model components has been critically verified and compared to analytical solutions in [25], two examples with particular relevance to PBFAM melt pool modeling, will be recapitulated in the following. The results of The only difference between the two variants is that the evaporative heat loss term (5) has been neglected in the second variant according to Figure 7 . In order to end up with a comparable effective energy input into the system and a comparable level of peak temperatures in the melt pool, the laser power had to be reduced by a factor of 200 for this second variant. On the one hand, this procedure demonstrates the importance of this heat loss term. On the other hand, it allows to study different physical phenomena and resulting melt pool characteristics, namely either melt droplets ejected from the melt pool ( Figure 7 ) or gas bubbles trapped in the melt pool ( Figure 6 ). First, these two scenarios allow to verify the robustness of the proposed SPH formulation in representing highly dynamic changes of the liquid-gas interface topology involving the coalescence and separation of interface segments. On the one hand, these are processtypical scenarios when spatter or gas inclusions are created, e.g., in the keyhole regime of PBFAM. On the other hand, an accurate and robust description of such interface phenomena is typically rather challenging for mesh-free discretization schemes such as SPH, which underlines the robustness of the present formulation. Second, these two variants also allow to study the different physical phenomena underlying the mechanisms of spatter generation respectively gas bubble inclusion, thus gaining detailed understanding of process dynamics and instabilities in PBFAM. Figure 6 displays the system behavior of variant 1. Accordingly, the interaction of surface tension and recoil pressure forces leads to oscillations of the liquid-gas interface with maximal amplitudes slightly above the bottom of the keyhole. Once, the amplitudes of these oscillations are large enough such that the opposite keyhole walls gets into touch a liquid bridge forms, which effectively encloses the gas material at the keyhole base -a gas inclusion is formed. Note again, that this example mainly aims at demonstrating the robustness of the proposed formulation in representing the highly dynamic surface tension-recoil pressure interaction and eventually the formation of a gas bubble. Of course, the employed phenomenological recoil pressure model does not explicitly resolve the high-velocity vapor jet that would arise from the keyhole base in the real physical system. This vapor jet might considerably influence the described creation mechanism of a liquid bridge, and potentially burst it before a closed gas bubble is created. In Section 3.3 first steps towards a high-fidelity melt pool model that explicitly resolves these effects of evaporation and gas dynamics will be presented. Also first steps towards the representation of mobile powder particles in the melt pool model have been made [80].

FIGURE 8
3D line melting problem with explicitly resolved powder particle geometry: resulting melt pool shape and final topology of solidified surface (post-processed) as well as temperature field in the range from 300 (blue) to 3700 (red) [25].
Eventually, Figure 8 displays two configurations of a 3D line scanning example with spatially resolved powder particle geometry. Accordingly, surface tension forces dominate volumetric forces such as gravity at the considered length scales and smoothen out the original particle contours almost immediately after melting. The peak temperatures at the melt pool center exceed the boiling temperature, and the resulting recoil pressure forces foster the creation of a depression at the melt pool center.

Alternative approach: High-fidelity melt pool model with resolved evaporation dynamics
In order to describe mechanisms of defect creation, e.g., creation of gas inclusions as discussed above, more accurately and to generally consider evaporation-induced gas/vapor flows and thereby induced material re-distribution dynamics (e.g., powder particle entrainment by gas flow), a model is required with spatially resolved gas and vapor phase as well as liquid-vapor phase transition. Our ongoing research work focuses on the development of such an alternative, high-fidelity melt pool model, which differs from the SPH model described above by the following main aspects: (i) evaporation, i.e., the phase transition from liquid to vapor/gas phase is explicitly resolved, and the associated (recoil) pressure jump and heat transfer across the liquid-gas interface follow consistently from the balances of mass, momentum and energy instead of employing the phenomenological models (4) and (5); (ii) an Eulerian discretization scheme (viz. the finite element method, FEM) is applied, which requires a tracking scheme (viz. the level set method) for the position of the (diffuse) liquid-gas interface; (iii) a truly incompressible instead of weakly incompressible flow is considered in the liquid and gas phase. For this approach, model equations (3) have to be supplemented by a transport equation for the level set method that additional allows for the liquid-vapor phase transition. Moreover, the continuity equation (3a) takes on the form of a constraint equation in the truly incompressible case. The pressure can be identified as associated Lagrange multiplier, and the densities and of the liquid and gas phase are constant and a priory known material parameters, i.e., the equation of state (3d) is not needed anymore in this case. In sum, the problem is described by a set of six equations, namely the incompressibility constraint equation, the Navier-Stokes momentum equation (3 components), the energy equation, and the level set transport equation, for the six primary unknowns given by velocity (three components), pressure , temperature and a level set function . The code implementation of the high-fidelity melt pool model 1 strongly relies on the software package adaflo [81] and the deal.II finite element library [82,83].   Figure 10 ) is compared to a variant of this FEM model, where the phenomenological recoil pressure model (4) is employed instead. For simplicity, and to isolate the effect of interest, the temperature field (relying on an approximate analytical solution) and the resulting melt pool shape are prescribed a priori for this example instead of solving the energy equation. Moreover, at the current state of the ongoing model development, the entire gaseous phase is assumed to consist of metal vapor, i.e., no local phase fractions of vapor and inert gas with correspondingly averaged material properties are considered. Please note also, that the parameters of the two variants have been chosen to lie in a comparable range, but due to the current state of model development a direct quantitative comparison of the results is not possible. Given these restrictions, the melt pool shapes in Figures 9 and 10 are in a reasonable qualitative agreement. For both variants, the evaporation-induced recoil pressure results in a deep depression at the center of the melt pool and melt spatter ejected from the pool at its lateral boundaries. Per underlying model assumptions, only the high-fidelity variant according to Figure 10 can predict the vapor velocity field resulting from the evaporation dynamics (visualized by white arrows).

Experimental validation
Again, drawing conclusions from the computational melt pool model requires experiment-based calibration, i.e. inverse identification, of the model parameters and validation of results. These critical steps will be performed as part of our future research work based on the expertise of the MIT team and their infrastructure in terms of process control and metrology instrumentation following two principle routes -in-situ optical process monitoring and ex-situ, post-process analyses.
Optical monitoring techniques of PBFAM have been successfully employed to detect defect creation and correlate process signatures to final part properties. For example, defects, spatter, melt track dimensions and even melt pool instabilities can been detected in-situ by identifying temperature developments in the laser scan path [84][85][86][87][88][89][90][91][92][93][94][95]. The team at MIT has developed a midwave infrared (MWIR) imaging technique, that overcomes limitations of previous techniques such as distortions and chromatic aberrations due to the usage of f -lenses or perspective distortion and limited depth of field due to the typical off-set viewing position [96][97][98]. The technique provides a high-resolution temperature field by capturing thermal emissions from the build surface. It applies a unique aperture division multiplexing (ADM) method, where the input aperture of a single large-diameter lens is sub-divided to bring the beam paths of the laser and the optical monitoring signal to focus on the build platform [99].
Calibration of model parameters can be conducted with e.g., simple point/line heating and melting experiments with varying laser powers. This can provide a good basis for model parameters such as laser absorptivity of the melt surface or temperaturedependence of the surface-tension, since the thereby induced Marangoni flow will considerably influence the temperature distribution in the melt pool. Also for the high laser power regime, where evaporation has considerable influence on the temperature profile and the generation of spatter, which is also identifiable via MWIR, point melting experiments are considered suitable to identify model parameters associated with the thermo-and fluid-dynamical phenomena of evaporation. Of course, process temperature signatures are also suitable for validation of actual powder melting simulations.
As second means for validation, ex-situ data will be considered. For post-process characterization, laser confocal microscopy of the sample surface, cross-sectional optical microscopy and scanning electron microscopy (SEM), as well volumetric imaging via micro-computed tomography (CT) can be employed. By this means, melt track morphology (e.g., shape, height, width), defects such as voids (size, shape, position, frequency), and surface topology can be used as metrics for model validation.

MACROSCALE THERMO-SOLID-MECHANICAL MODELING
The TUM authors recently proposed a thermo-mechanical finite element model aiming at the prediction of residual stresses and thermal distortion in the partscale simulation of metal PBFAM processes [13,100]. The most important model equations as well as exemplary results are recapitulated in the following.
An alternative thermo-mechanical finite element model for PBFAM partscale simulations has been proposed by N. Hodge and the team at LLNL [5,101]. In particular, the PBFAM example in the microstructure studies of Section 5.2 has been simulated using this model and its code implementation in the in-house parallel finite element code Diablo [102] at LLNL. Recent extensions of the underlying numerical schemes focus on the critical aspects of computational efficiency and accuracy by adequately addressing the spatial and temporal multiscale nature of these processes [12,103]. For further details on this alternative partscale modeling approach the interested reader is referred to the aforementioned references.

Model equations
The thermo-mechanical problem consists of the dynamic heat equation coupled with the static balance of linear momentum: with the primary variables temperature and displacement . The magnitudes of strains and rotations arising from typical PBFAM processing conditions can be assumed as small, hence commonly the theory of linear continuum mechanics in combination with the engineering strain tensor as metric for deformation are employed. The structural material law relating those with the Cauchy stresses = ( ( ), ) will be detailed Section 4.1.2. The heat flux is specified by Fourier's law of heat conduction, The material parameters appearing in the thermal problem, namely volumetric heat capacity and heat conductivity , in general depend on temperature and phase as discussed in the following Section 4.1.1. Finally,̂ represents a volumetric laser beam heat source, which is modeled according to [2].

Temperature-and phase-dependent parameters
To model the three relevant phases powder, melt and solid, in a first step the liquid phase fraction is defined according to where and represent the liquidus and solidus temperature. In a next step, the irreversibility of the powder-to-melt transition during melting is captured via the consolidated fraction In a last step, the final phase fractions of the powder ( ), melt ( ) and solid ( ) material are computed following the relations: = , Eventually, these phase fractions can be used to interpolate an arbitrary material parameter according to the following scheme: where interp is the interpolated parameter value and , and are the single phase parameters. This approach is applied to the thermal conductivity and heat capacity . The phase-and temperature-dependent formulation of mechanical material properties is presented in the next section.

Mechanical constitutive law
An iso-strain (Voigt-type) homogenization, assuming equal strains in all three phases, is applied in the following. Accordingly, the stress of the mixture is given by a weighted sum of the individual contributions, a procedure that is in fact similar to the interpolation scheme (15): In general, the strain (8) in a single phase ∈ { , , } is calculated from an additive decomposition of the strain tensor although not all terms will be utilized for each phase. The total or kinematic strain on the left-hand side of (17) is purely displacement-dependent and for all phases given by the kinematic relation (8). The first term on the right-hand side of (17) is the elastic strain , , which leads to a stress in each phase when considering a linear hyper-elastic material with fourth-order elastic constitutive tensor C , e.g., according to Hooke's law with Young's modulus and Poisson's ratio (assumed to be equal in all phases) as independent constitutive parameters. The plastic strains , are only relevant for the solid phase, and can be calculated using standard approaches, e.g., an incremental problem formulation in combination with a return mapping algorithm. For simplicity, this contribution will not be considered in the examples discussed below. The strains , due to thermal expansion are assumed equal in all phases and read where is the coefficient of thermal expansion. Critically, an inelastic reference strain ref, , which is only relevant for the solid phase, is proposed in rate form according to: Note, that the reference strains only change when the solid phase fraction increases according tȯ > 0 (first case in (20)), i.e., for temperatures ∈ [ ; ] in the phase change interval and negative temperature rateṡ < 0. Note also that the second case in (20)  An elastic constitutive law with low stiffness values (i.e., , , ≪ ) as applied to powder and melt leads to small stresses yet considerable total strains in these phases. In this context, the reference strains according to (20) ensure that these strains do not translate into stresses during solidification. For the special case that kinematic , plastic and thermal strains are constant during solidification, which approximately holds if the phase change interval − is sufficiently small, it can easily be verified from (17) and (20) that the elastic strain, and thus the resulting stresses, in the evolving solid phase vanish.

Remark:
One of the main assumptions underlying the present and most existing thermo-mechanical PBFAM models is that mechanical stresses in the (open-surface) powder and melt phase domains are negligible. This behavior is approximated by applying a simple elastic constitutive law to these phases, with stiffness parameters that are considerably lower as compared to the solid phase, i.e., , ≪ . In practice, this approximation turns out to result in moderate, i.e., limited, strains, since the deformation of these powder and melt domains is mostly kinematically controlled by the motion of the significantly stiffer solid phase domains, thus yielding only small stress contributions as desired. Moreover, as compared to approaches exactly satisfying the zero-stress assumption in powder and melt, no additional means (e.g., extended finite element method, immersed boundary method, etc.) are required for representation of discontinuities (e.g., jumps in stresses) inside elements or for mesh movement in the geometrically nonlinear case. Note, the assumption that thermal strains exist also in the powder and melt phase, and are equal to thermal strains in the solid phase, has been made for simplicity here. This assumption is neither necessary nor has it a significant influence on the resulting residual stresses due to the low stiffness of these phases and the definition of reference strains (20), which ensure that newly created solid material is stress-free.

Exemplary simulation results
In [13], the accuracy of the proposed thermo-mechanical model and the influence of the individual model components has been critically verified by means of elementary test cases, partly with analytical solutions.
In the following, two examples with direct relevance for macroscale modeling of PBFAM processes will be briefly recapitulated. The first example consists of a solid base plate (confined by horizontal solid line) and ten melt tracks (separated by dashed horizontal lines) successively deposited on top of each other. Here, the length of the (horizontally centered) melt tracks has been chosen as half of the length of the base plate, and the laser is scanning from left to right in each track. The resulting normal stresses in x-direction, i.e., in horizontal direction, for different snapshots in time are illustrated in Figure 11 . As characteristic feature of the stress distribution, a horizontal band pattern can be identified, with stresses alternating between positive (tensile) stress values in the upper part and stress values close to zero in the lower part of a each melt track. In the following, the creation of this track-wise band pattern will be explored exemplarily during processing of track 7. The top left snapshot of Figure 11 represents a configuration, where track 6 has been finished and cooled down already. The subsequent snapshots (from left to right) in the first and second row of Figure 11 show the deposition of track 7. To visualize the melt pool size, the temperature isoline corresponding to the solidus temperature is depicted in these snapshots (solid black line). While no stresses occur in the powder material of layer 7 in front of the laser, the thermally induced volume expansion during heating leads to negative (compressive) stresses in the solid material of the previously processes track 6, mostly pronounced in the direct vicinity of the

FIGURE 12
Processing of a solid and hollow pyramid. Discretization with non-matching meshes via mortar mesh-tying [13].
-isoline. As desired, these stresses rapidly drop to zero in the narrow temperature interval ∈ [ ; ] such that no visible stresses remain in the melt pool domain. This strong gradient between vanishing stresses in the melt pool and high compressive stresses in the solid material beneath remains after solidification and is superimposed by additional tensile stress contributions due to the thermally induced volume shrinkage during cooling. After cool down (see snapshot at bottom right of Figure 11 ), this process results in high tensile stresses in the upper, re-molten part of track 6, and stresses close to zero in its lower part. The same characteristics are observed in the previously processed tracks below. Even though the base plate has the same stiffness and similar support conditions (i.e., spatially fixed only at its bottom surface) as the solidified tracks above, this characteristic temperature gradient is much more pronounced for the first track, i.e., the highest overall tensile stresses occur in the first track, accompanied by compressive stresses of comparable magnitude in the base plate below. This observation can be explained by the fact that the base plate is initially stress-free, while solidified melt tracks are subject to tensile stresses after cooling down (e.g., snapshot at bottom right of Figure 11 ), which partly compensate the compressive stresses arising from thermal expansion when processing the subsequent track above.
As important feature of the proposed macroscale modeling framework, successive layers of material are connected by powerful mortar mesh-tying schemes [104,105]. These schemes allow for non-matching meshes between these layers while preserving optimal convergence rates in the 2 -norm, which is typically not guaranteed by alternative, e.g., collocation-type, mesh-tying approaches. In particular, the combination of complex geometries, e.g., tapered parts with strongly changing cross-section dimensions, and the need for layer-wise mesh definition makes the discretization of PBFAM parts challenging when using standard discretization schemes. In Figure 12 , a geometry of this type is illustrated, consisting of a solid pyramid with curved edges and an upside-down hollow pyramid. Note that the meshes used for discretization of the pyramids are non-matching between successive layers, and also non-matching between the pyramids and the base plate. The flexibility arising from this non-matching discretization approach allows to have very regular meshes with undistorted and almost equally-sized finite elements, which does not only decrease the effort and time for mesh generation but also the number of finite elements, i.e., degrees of freedom, to achieve a certain approximation quality, i.e., to lie below an acceptable discretization error level.
Again, the visualization of von Mises stresses in Figure 12 (a) reveals that highest residual stresses occur in the first layer above the base plate. The displacement magnitudes associated with the coupled thermo-mechanical problem are illustrated in Figure 12 (b). Accordingly, the thin walls of the hollow pyramid exhibit a strong thermal distortion during cool down.

Experimental validation
Our future research work will focus on experimental validation of the macroscale model. Quantification of residual stress and distortion, respectively, require measurements of the lattice distortion of solidified material and the three-dimensional component shape. Both are intrinsically challenging, and therefore much research has relied on fabrication and metrology of geometric test artifacts. As such, early work by Kruth and colleagues studied the influence of scan pattern (e.g., unidirectional vs. raster scan) on residual stress development in small arch-shaped geometries [106]. The coupled influence of support structure placement and geometry on deformation and residual stress has also been studied. In situ measurements of stress buildup in LPBF were performed using a build plate instrumented with strain sensors [107], or, alternatively, by embedding strain-sensing fiber optics in SLM components [108]. High energy X-ray and neutron diffraction scattering enable local, volumetric probing of residual stress, and can be employed ex situ (e.g., after printing) as well as during heat treatment. The sample geometry must ensure beam penetration; neutron diffraction is attractive due to the higher penetration depth through metallic materials [109], where X-ray techniques are preferred for thin film geometries. Last, three-dimensional imaging of surface and shape, such as by digital image correlation (e.g., [110]) and computed tomography can be compared to simulations of component deformation.
Foundational work by N. Hodge and the team at LLNLHodge [111] demonstrates tight agreement between a continuum model and residual stress measurements by digital image correlation and neutron diffraction. Techniques for stress mitigation and active control are especially necessary when fabricating components where geometric accuracy of as-printed features or predictable stress-strain behavior is paramount, both in LPBF and directed energy deposition (DED) methods.

MICROSTRUCTURE MODELING
The authors recently proposed a microstructure model for phase fraction prediction in PBFAM of Ti-6Al-4V [112]. The most important model equations and exemplary results are recapitulated in the following.

Model equations
Instead of spatially resolving individual crystal or grain geometries, the model describes solid state transformations in Ti-6Al-4V via spatially homogenized phase fractions of the most relevant metallurgical specifies, namely the -phase, the stable s −phase as well as the metastable martensitic m -phase. As compared to spatially resolved approaches, this procedure offers the general suitability for partscale PBFAM simulations. The temperature field, which is required as input for the microstructure model, can e.g., be provided by macroscale thermal or thermo-mechanical PBFAM models according to Section 4. The composition of the solid phases , s and m are described by phase fractions X ∈ [0; 1] fulfilling for fully solidified material. Here, the stable and martensitic phases, s respectively m , are combined to the total alpha phase fraction X . For sufficiently slow cooling, only the stable phases s and will arise. The temperature-dependent equilibrium phase fraction X eq ( ), towards which the s -phase tends in the extreme case of very slow cooling rates |̇ | ≪ |̇ m ,min |, is described as Koistinen-Marburger law [113] according to: X eq ( ) = 0.9 for < ,end , Thus, for temperatures above the alpha-transus end temperature ,sta = 935 pure material, and for temperatures below the alpha-transus end temperature ,end = 848 a phase composition with 10% -and 90% stable -phase are expected under thermodynamic equilibrium conditions. For the second extreme case of very fast cooling rates |̇ | ≥ |̇ m ,min |, at which the diffusion-driven formation of X s is suppressed, a metastable Martensite pseudo equilibrium X eq m ,0 ( ) is stated [56,113,114]: Here, |̇ m ,min | is denoted as critical cooling rate, above which pure Martensite formation is observed. Note that this critical rate is not prescribed in our model but emerges naturally from the system dynamics. Note also that, according to this model, Martensite formation can only be initiated for temperatures below the Martensite start temperature m ,sta , and a maximal Martensite phase fraction of 90% can be achieved at room temperature ∞ = 293 . Finally, in the most general case of cooling rates that are too fast to complete the diffusion-driven s -formation before reaching the Martensite start temperature m ,sta but still below the critical rate |̇ m ,min |, i.e., a certain amount of stable s -phase has already been formed at m ,sta , a Martensite phase fraction below 90% is expected at room temperature. For this case, we propose to replace (23) by an effective pseudo equilibrium phase fraction X eq m ( ) accounting for the reduced amount of transformable -phase at presence of a pre-existing phase fraction X s : Eventually, the formation and dissolution of the three species is described in rate form, i.e., evolution equations are proposed with the following contributions to the total ratesẊ s ,Ẋ m andẊ : Here, e.g.,Ẋ → s represents the formation rate of s out of whileẊ s → represents the dissolution rate of s to . Since thephase fraction can be calculated from the continuity constraint (21), the first two rate equations in (25) are sufficient to predict the evolution of X s and X m . Based on the physical mechanism underlying the phase transformation, the formation and dissolution rates in (25) are classified as instantaneous vs. diffusion-based transformations, which are initiated through deviations of the current martensitic and total -phase fraction X m respectively X = X s + X m from the corresponding equilibrium values according to (24) and (22). While instantaneous transformations are modeled via constraint equations enforcing that phase fractions follow the associated equilibrium value, the diffusion-based creation of a phase out of a phase is modeled as a modified logistic differential equation [115] of the following form: Here, the factor ( − ) represents the driving force of the diffusion process in terms of transformable -phase and has a decelerating effect on the transformation during the ongoing diffusion process. The factor with leads to a transformation rate that increases with increasing amount of created s -phase, i.e., it has an accelerating effect on the transformation rate during the ongoing diffusion process. The factor ( ) represents the temperature-dependent diffusion rate of this thermally activated process, taking into account the temperature-dependent mobility of the diffusing species. Eventually, the choice of the exponent dependents on the specific physical mechanisms underlying the diffusion process. Specifically, the ratesẊ → m for Martensite formation out of andẊ m → for Martensite disolution into are modeled as instantaneous transformations, the remaining rates in (25a) and (25b) as diffusion-based processes according to (26). Further details are given in [112].

Exemplary simulation results
In [112], time temperature transformation (TTT) experiments are used for model calibration, i.e., to inversely determine the unknown model parameters (e.g., in the diffusion equation (26)) such that the deviation between phase fraction contour-lines from experiments and corresponding simulations are minimized in a least-square sense. In TTT experiments, the material is first equilibriated at high temperatures such that only the high-temperature phase, here the -phase, is present. Subsequently, the material is rapidly cooled down to a target temperature at which it is then held constant over time so that the isothermal phase transformation at this temperature can be recorded. Rapid cooling refers here to a cooling rate that is fast enough so that diffusion-based transformations during the cooling itself can be neglected and can subsequently be studied under isothermal conditions at the chosen target temperature. The procedure is repeated for successively reduced target temperatures. The emerging diagram of phase contour-lines over the × log( ) space is commonly referred to as TTT-diagram, which has been plotted for experimental measurements and model predictions of X s and X m in Figure 13 . Note that the depicted phase fraction values are normalized to the corresponding equilibrium value. In addition to the TTT experiments used for model calibration, continuous-cooling transformation (CCT) experiments have been considered in [112] for model verification. In CCT experiments, the microstructural probe is again equilibriated at high temperatures where only the high-temperature -phase is present. Afterwards, the probe is cooled down to room temperature ∞ = 293.15 at different cooling rateṡ CCT following precisely defined, time-continuous temperature profiles. Eventually, the evolving microstructure resulting from this cooling procedure is recorded on the × -space. It was demonstrated that the model prediction of the critical cooling rate |̇ m ,min | ≈ 410 ∕ above which pure Martensite formation was observed, is in very good agreement to experimentally measured values. In many existing models this critical cooling rate is enforced as ad-hoc criterion for Martensite formation. In our model this characteristic cooling rate is not prescribed but is a direct consequence of the energy and mobility competition between the microstructural species.
As first PBFAM-specific example, Figure 14 shows the simulated microstructure composition in the vertical center plane of a cube fabricated by selective laser melting. For this specific problem, the required thermal field has been provided by a macroscale PBFAM model developed and implemented at the LLNL [5,101]. To accurately predict the microstructure evolution in PBFAM, it is necessary to spatially resolve the correct laser path and layer/track dimensions instead of applying layer/track up-scaling or agglomeration approaches [3,5,9]. Since the partscale simulation of PBFAM with resolved laser path is still an open research question, the following example will be limited to a cube of side length 1mm consisting of 34 layers. According to Figure 14 , no stable -phase is formed due to the fast cooling rates in PBFAM. Instead, the -phase directly transforms into martensitic m -phase once the temperature falls below the Martensite start temperature m ,sta .

FIGURE 13
Simulation of TTT-diagram for the s -and m -phase, along with experimental data by Malinov [116] and Kelly [117]. Left: Contour-lines for X s ; Right: Contour-lines for X m . Contour lines are shown for the 1%, 5%, 45%, 55%, 95% and 99% normalized phase fractions. Three temperatures are marked in red and discussed in the analysis [112].

FIGURE 14
Simulated microstructure composition within the vertical center plance of a one millimeter-sided cube fabricated with SLM. The gray area above marks the layers that are not yet processed by the laser at the considered time. In the right figure, which depicts the -phase fraction X , the melt-pool is indirectly visible (the liquid phase fraction is not depicted) through the decreased -phase fraction in the right corner of the up-most layers due to the laser that has just previously scanned this plane in x-direction from left to right. In a similar fashion, the decreased amount X in the upper left corner of the right Figure stems from the heat of the laser that is already melting the subsequent track at the considered time [112].

FIGURE 15
Simulated microstructure distribution of s -and m -phase in the vertical center plane of a ten centimeter-sided cube resulting from rapid cooling via thermo-mechanical contact at the bottom and free convection at the remaining surfaces [112].
In PBFAM parts of practically relevant size, the repeated thermal cycles in the process lead to temperature evolutions where material points in the center region of the part remain at elevated temperatures for considerable time spans such that also the creation of -phase, either directly out of the -phase or via diffusion-based dissolution of Martensite at elevated temperatures (both effects are captured by the proposed model) can be expected. To investigate the influence of the increased thermal mass of larger parts on the centimeter scale, a final example was considered in [112], where the quenching process of a cube with 100 mm side length was studied. The cube with a homogeneous initial temperature of 0 = 1300 is in thermo-mechanical contact with a cold base-plate at its bottom surface (fixed temperature of bp = 300 at the bottom; effective heat-transfer coefficient of tc = 5 ⋅ 10 5 2 at the (upper) contact interface with the cube) and subject to a convective boundary conditions (heat-transfer coefficient c = 1000 2 ) on its remaining (free) surfaces. The microstructure evolution resulting from the rapid cooling induced by these boundary conditions is illustrated in Figure 15 .
Accordingly, a several millimeter thick Martensite coating at the surfaces is observed, which is especially pronounced at the corners of the cube that are characterized by higher cooling rates and lower temperature levels. These qualitative findings are in good agreement with experimental results for quenching of Ti-6Al-4V as well as additively manufactured parts [118,119]. Due to the higher heat-transfer coefficient, a thicker Martensite coating is found at the bottom surface of the cube as compared to the remaining free surface areas. As expected, the core of the cube, where cooling rates are lower and the temperature remains at elevated levels for longer times, is composed of stable s -phase. Moreover, in [112] it was demonstrated that heating of the base plate, either in the quenching or in the PBFAM example, is an effective means to suppress local Martensite formation.

Experimental validation
In addition to the experimental model validation via TTT-and CCT-diagrams as presented above, further steps of AM-related experimental microstructure characterization will be considered in our future research work. Well-established techniques for metallurgical characterization are typically used to identify the microstructural details of components made by LPBF and other metal AM processes. Optical microscopy of surfaces, typically etched to enhance the appearance of grain boundaries, is a routine method for identifying grain size and shape. Directional reflectance microscopy, combined with image processing, has been used for automated size segmentation [120]. Electron microscopy based methods such as electron backscatter diffraction (EBSD) are most commonly used to assess grain structure and orientation, again by surface analysis. Careful preparation of the sample surface for imaging, often by mechanical sectioning and polishing, is necessary and is commonplace in laboratories. Higher-resolution information, down to the nanoscale and atomic scale distribution of precipitates and alloying elements, can be obtained by transmission electron microscopy, atom probe tomography, and other techniques [121]. These experimental methods are used to understand the relationships between process parameters, microstructure, and ultimate performance. For example, Kohnen and coworkers demonstrate how scan velocity influences EBSD-measured grain size and texture, and in turn the performance of tensile test artifacts fabricated via SLM [122] (also see, e.g., [123]).

CONCLUSION AND OUTLOOK ON OVERALL MODELING APPROACH
In this article, an overview was given on the authors' recent developments in the modeling of powder bed fusion additive manufacturing (PBFAM) processes on different length scales, specifically with respect to mesoscale powder modeling, mesoscale melt pool modeling, macroscale thermo-solid-mechanical modeling and microstructure modeling. Beyond the insights gained from these individual models, also their interplay in an integral modeling approach will be a central aspect of future research.
The long-term objective of the presented mesoscale modeling approaches is to model the entire process chain from powder to part, i.e., the powder spreading as well as the subsequent laser melting and solidification process in PBFAM. First, the powder modeling approach according to Section 2 aims at the identification of key feedstock properties allowing to characterize the spreadability of new powder materials on basis of simple powder-rheological experiments. Second, the powder spreading simulations are intended to support the development of improved powder spreading strategies, e.g., novel spreading tools and kinematics, to enable controlled and repeatable spreading conditions that allow for high layer quality in terms of surface uniformity as well as high and spatially constant packing density. Specifically, a focus will lie on fine-grained, highly cohesive powders, for which a controlled spreading is challenging due to the dominating role of cohesive forces and thereby reduced flowability, but which are desirable for processes such as electron beam melting (EBM), selective laser melting (SLM) or binder jetting (BJ) for reasons of increased geometrical resolution, lower material costs and, in case of binder jetting, reduced sinter times. Taking the layers provided by the powder spreading model as input, the melt pool simulations of Section 3 aim at correlating these layer properties with melt pool dynamics respectively instabilities, and finally with resulting defect characteristics, e.g., size, shape and frequency of pores. Specifically, the question shall be answered how melt pool stability can be fostered by improved opto-thermo-mechanical characteristics of high quality powder layers.
On the macroscale, thermo-mechanical models according to Section 4 are intended to predict temperature fields, residual stresses and thermal distortion on the scale of design parts. The predicted temperature fields can be taken as input to predict metallurgical phase fraction evolutions using (homogenized) microstructure models as presented in Section 5. The knowledge of these phase fractions can be exploited to formulate microstructure-informed, macroscale constitutive laws, which are calculated by spatially homogenizing the material properties of the individual species, and which allow for more accurate residual stress predictions when integrated into thermo-mechanical PBFAM models. Due to the extreme spatial and temporal temperature gradients during PBFAM, which fosters the formation of non-equilibrium (e.g., martensitic) phases, and the considerably different constitutive behavior of these diverse phases (e.g., ductility of -phase vs. martensitic -phase in Ti-6Al-4V), this approach is expected to significantly increase model accuracy.
Also an information transfer between macro-and mesoscale can be beneficial for overall modeling. On the one hand, macroscale models can provide improved thermal boundary conditions for representative volumes at different locations of the design part, which are considered for mesoscale melt pool simulations. On the other hand, the mesoscale powder and melt pool models may be exploited to derive improved effective continuum properties (e.g., effective opto-thermal properties of powder phase, anisotropic thermal conductivity as model for convective heat transfer in the melt, etc.) for the macroscale model.
Many of the parameters underlying the aforementioned models cannot be taken from standard databases due to the extreme processing conditions (e.g., extreme spatial and temporal temperature gradients) in PBFAM and the high sensitivity/variability of certain model parameters with respect to physical state variables (e.g., temperature, pressure, etc.), environmental conditions and material imperfections. In our ongoing research, high-resolution measurement data from powder spreading and laser melting experiments as described in Sections 2.3 and 3.4 will be combined with probabilistic schemes for inverse parameter identification and uncertainty quantification to identify the unknown model parameters under process-relevant conditions. In addition to parameter identification/calibration, the underlying methods for inverse analysis will also be used to directly (i.e., without repeated and computationally expensive forward simulations in a 'trial-and-error' manner) answer practically relevant inverse questions, i.e., the question which input is required to achieve a desirable output. An example in this context is the geometrical compensation of thermal distortion by solving the inverse design task of finding an optimal initial geometry Given the multiscale nature of metal PBFAM processes and the complexity arising from competing physical mechanismen and various types of defects on these different length scales, an integrated modeling and simulation approach unifying information from macro-, meso-and microscale models has the potential for virtual process optimization and part qualification, thus promoting broader industrial adoption. model equations, the code implementation and conduction of numerical simulations w.r.t. the powder model. M. Schreter contributed to the conceptualization, derivation of model equations, the code implementation and conduction of numerical simulations w.r.t. the (FEM) melt pool model. N. Hodge contributed to the conceptualization, derivation of model equations, the code implementation and conduction of numerical simulations w.r.t. the microstructure model. A.J. Hart contributed to the conceptualization and design of the different experimental approaches, and supervised work at MIT. W.A. Wall contributed to the conceptualization of the different modeling approaches, and supervised work at TUM.