Non-isothermal Phase-Field Modeling of Heat–Melt–Microstructure-Coupled Processes During Powder Bed Fusion

Modeling and simulation of powder bed fusion (PBF) remain a great challenge due to the sophisticated and interactive nature of underlying physics. A unified scenario considering interactions among the heat transfer, melt flow dynamics and microstructure evolution (noted as “heat–melt–microstructure-coupled processes”) is therefore essential for a thermodynamically consistent description and thus reliable microstructure prediction. In contrast to the state of the art, where either individual aspects are considered or the thermal history is taken as input from separate numerical scheme, we propose in this work a unified non-isothermal phase-field model for the heat–melt–microstructure-coupled processes during PBF. Simulations on a stainless steel 316L powder bed demonstrate that the model can reproduce well-observed features, but also help to discover new in-process phenomena and reveal the mechanism of the defect formation. Based on massive simulation results, we also present the densification map with respect to beam power and scan speed, and have classified the regions of the parameter combination by the distinct resultant morphology.


INTRODUCTION
Distinguished from conventional subtractive methodologies, additive manufacturing (AM) implies layer-by-layer shaping and consolidation of feedstock to arbitrary configurations. Powder bed fusion (PBF) is a family of the common techniques that actualize the AM fabrication in a fashion of sequentially layer-by-layer powder spreading and subsequent binding driven by beam scans, such as selective laser melting/sintering (SLM/SLS) and selective electron beam melting. [1][2][3] Although PBF is conceptually simple, the processes are not yet fully controllable, reproducible, or predictable. 1 As a result, its wide application is currently limited mostly by quality issues of the end-use products, e.g., insufficient hardness and fracture resistance due to the deformation, cracks, and porosity. In recent years, experimental efforts focusing on optimization of the process parameters, mostly related to beam power, scan speed, and powder materials, have achieved significant progress. Nevertheless, these investigations often suffer from the trial-anderror principle, which is time-consuming and not sufficient to provide transferable fundamental insights into the process-microstructure-property relationships.
To this end, modeling and simulation of PBF on the mesoscopic scale of powders indubitably play a crucial role in successful manufacturing, since it allows the direct linkage between the process parameters and the resultant microstructure that determines the macroscopic performance of the final product. 1 However, the underlying physics are complicated and cover a broad range of time and length scales. These physics include but are not limited to beam-powder bed interaction, heat and mass transfer, phase transitions, melt pool behaviors, and structural relaxation. 4 To date, various methods have been applied or extended for the simulation of PBF targeting at specific issues. Computational fluid dynamics (CFD), implemented with finite element method (FEM) and finite volume method (FVM), is a common approach to model the melting process as well as melt flow dynamics. Among existent CFD models, Navier-Stokes equations of incompressible Newtonian (INS) are widely utilized to describe the behavior of melt flow, which is usually associated with heat transfer, inertia force, capillary effect (also known as wetting) and/or thermocapillary effect (also known as Marangoni flow), and vapor pressure. The volume-of-fluid approach is additionally implemented to identify the evolution of the liquid-gas interface. Khairallah et al. presented the first three-dimensional (3D) model of a SLM single scan using the CFD code implemented with a hybrid of FEM and FVM. Arbitrary Lagrangian-Eulerian has also been implemented to deal with the mesh vertices which may flow with the fluid (Lagrangian) or be fixed (Eularian). [5][6][7] A well-resolved mesoscopic powder bed geometry was also created in order to reproduce the experimental setup with least approximations. This model was later extended with recoil pressure and Marangoni convection, 8 and denudation effects of metal powders 9 for laser-based PBF. Körner et al. reported the first numerical lattice Boltzmann model (LBM) to analyze melting processes on the powder scale during SLM. 10 This model was extended further with physics, such as the phase transitions and capillary forces .11,12 Besides the efforts on melt pool behaviors and pore structure evolution, generation of grain structure is also a highly concerned topic in modeling and simulation of PBF, due to its significant influence on the mechanical and functional (e.g., magnetic, electric or optical) properties of the final product. Commonly used methods include cellular automata (CA). [13][14][15][16] Rai et al. reported a CA-LBM coupled model, which simulates two-dimensional (2D) re-solidification under a uniform thermal field. 13 Recent publications have shown that the CA approach can also be successfully employed for the 3D simulation of coupled re-solidification. [15][16][17] Nevertheless, the above-mentioned works are subject to strong simplifications and segregated modeling schemes, considering either only selected aspects or with the thermal history taken as input from separate numerical schemes. In fact, the heat transfer and thus also the melt pool is sensitive to the powder bed morphology, and has a strong local impact on melting, re-solidification, pore/bubble evolution and grain growth. A unified scenario considering interactions of the mentioned processes is thereby essential. Only in this way can one arrive at a thermodynamically consistent modeling of the PBF system and thus reliable microstructure prediction. One promising approach is phase-field modeling. [18][19][20][21] Phase-field models apply order parameters (OP) to represent the microstructure (pore/liquid/solid phases, grains, or even the solidstate phases), and can involve natural thermodynamic quantities like temperature (or internal energy) and mechanical fields. The OPs take different constant values in different regions and have rapid spatial variation across the interfaces. The microstructure and its evolution are thus recaptured by the distribution of the OPs in space and time, with no need for interface tracking. Thanks to the features such as generality of the theoretical framework, straightforward numerical implementation and coupling of multi-physics, phase-field modeling have been extensively studied for solidification, 19,21 grain growth, [22][23][24] sintering, 25 and fluid dynamics. [26][27][28] Even though these models are mostly proposed for isothermal cases with specific formulation of free energies, they provide a solid and rich background for phase-field modeling of non-isothermal cases and for application in PBF modeling.
Surprisingly, very few phase-field efforts have been made for the AM process, and they are mostly limited to individual microstructure aspects. Taking the temperature gradient and solidification velocity as parameters in established isothermal phase-field models, Gong and Chou 29 and Sahoo and Chou 30 simulated the variation of the dendrite morphology in a single scan. Nevertheless, there is still debate as to whether and to what extent a dendritic structure forms, since PBF-AM involves a very high solidification rate and solidification with a planar front may occur. Krivilyov et al. 31 presented phasefield simulations based on their multiphase flow model to study powder consolidation behavior, but with no explicit thermal information of AM. Zhang et al. 32 directly adopted the phase-field model constructed under an isothermal condition for SLS, while the heat transfer and the coupled nonisothermal microstructure evolution were not addressed. Yang et al. 33 proposed a thermodynamically consistent non-isothermal phase-field model for SLS, and proceeded with 3D simulations with a mesoscopic powder bed geometry. This model recapitulates processes during SLS, such as a temperature field with a high gradient, mass transfer through partial melting, and diffusion, particle/grain necking and coarsening. For the simulation of PBF, melting is dominant and the melt flow has a large influence on the pore and surface morphology. In 2018, Lu et al. 34 made the first attempt to construct a unified phase-field simulation of melting, re-solidification and grain structure formation under an evolving temperature field in the same spatially resolved powder bed. Their 2D simulations demonstrated that the model can reproduce many important microstructure phenomena observed experimentally. However, the model disregards the fluid dynamics of the melt pool.
In recent decades, the phase-field method has been successfully applied to solving fluid-phase-coupled problems involving a transient interface and/or topological variations (e.g., fluid jet, droplets and large-deformation waves) by the so-called Navier-Stokes-Cahn-Hilliard (NSCH) system 35,36 and the Navier-Stokes-Allen-Cahn (NSAC) system. 37,38 This shows the potential to model and simulate the PBF system in a unified and thermodynamically consistent way, which, however, has not been incorporated in the existent phase-field works.
In this work, we propose a unified non-isothermal phase-field model considering coupled processes among heat transfer, melt flow dynamics, and microstructure evolution (noted as ''heat-melt-microstructure-coupled processes''). The model is constructed from the entropy level following our former work 33,39 explicitly considering the contributions from heat transfer, configurational and gradient terms, which results in a free energy density functional receiving the influences from pore/liquid/solid phases as well as grains inside the solid phase, and the local temperature field. A coupled kinetic system for OPs (including the pore-substance, solid-liquid inside substance, and grains inside solid), melt flow dynamics and the thermal evolution under non-isothermal conditions are also derived. As examples, we perform 2D simulations of the PBF single scans on a stainless steel 316L (SS316L) powder bed. The influences from the beam power and scan speed on key features, such as the temperature field, melt pool geometry and interior process including grain growth and pore morphology evolution, are also discussed.

NON-ISOTHERMAL PHASE-FIELD MODEL
In this model, a conserved OP q is employed to represent the substance and atmosphere/pores; a set of non-conserved OPs / s and / l are employed to represent the solid and liquid phase, respectively; and a series of non-conserved OPs fg j g are employed to represent the orientation distribution among the powders/grains. According to this scenario ( Fig. 1), it is anticipated that none of the non-conserved OPs are valued when q ¼ 0, and when q ¼ 1, only one of / s and / l equals unity. This becomes the first constraint of OPs to restrict the existence of liquid as well as solid within the substance: ð1Þ On the other hand, fg j g have zero value when q ¼ 1, / l ¼ 1 and / s ¼ 0, indicating that the substance is in the liquid phase and thereby no grains should be present. When it is inside the solid substance (q ¼ 1, / l ¼ 0 and / s ¼ 1), for a specific grain j, only g j takes the unit value while the others g j Ã ; j Ã 6 ¼ j are equal to zero. This becomes the second constraint of OPs to have polycrystal inside the solid substance: Following this scenario, the entropy density functional S is constructed based on our latest work in Refs. 33 and 39: U ht ðq; f/ i g; fg j gÞs ht ðe ht Þ þ s cf ðq; f/ i g; fg j gÞ Â À s grad ðrq; fr/ i g; frg j gÞ Ã dX: ð3Þ with: s grad ðrq; fr/ i g; frg j gÞ ¼ where the entropy density s includes the contributions from thermal (s ht ), configurational (s cf ), and gradient (s grad ) contributions. j q , j / and j g represent the gradient constant. Here, de ht represents the heat change within X. When there is no volume variation, i.e. X is a constant, we have the internal energy density change de ¼ de ht according to the first law of the thermodynamics. The monotonic interpolating function U ht maps the heat contributions to the regions with a certain value of OPs q, f/ i g ði ¼ l, sÞ and fg j g ðj ¼ 1; 2; . . . ; NÞ. Correspondingly, the internal energy density functional E of a finite subdomain X within the system can be formulated as: ð4Þ where the internal energy density e consists of heat (e ht ) and potential (e pt ) terms. To relate the entropy density functional to the free energy density functional F, the Legendre transformation is performed 40 : inf e eðq; f/ i g; fg j gÞ À Tsðe; q; f/ i g; fg j gÞ Â Ã dX ¼ Z X U ht ðq; f/ i g; fg j gÞf ht ðTÞ þ f loc ðT; q; f/ i g; fg j gÞ Â þ f grad ðT; rq; fr/ i g; frg j gÞ Ã dX: The term f ht ¼ inf e ht ½e ht À Ts ht ðe ht Þ represents the contribution from the transformed heat term e ht , which leads to the relationship dðf ht =TÞ ¼ e ht dð1=TÞ: 41 Here, the heat term is calculated by the temperature change and the latent heat from the solid-liquid phase transition, i.e., e ht ¼ U ss c r ðT À T M Þ þ U l L; 33,39,41 where c r ¼ U s ðc s À c at Þ þ U l ðc l À c at Þ is the relative specific heat, setting the internal energy of the atmosphere/pores at a reference temperature T M of the system as the zero potential. c l , c s and c at are the volumetric specific heat of the solid, liquid and the atmosphere/pores, while U ss , U l and U s are the interpolating functions to indicate the substance, liquid and solid phase, respectively. L represents the latent heat of the liquid-solid transition. In this fashion, f ht can be formulated as: ð6Þ On the other hand, we simply formulate the local term f loc ¼ e pt À Ts cf into a polynomial, presenting the multi-well-shape with barrier heights linearly correlated with temperature, by adopting the same form for both e pt and s cf . As for a general PBF system, we specifically use three polynomials to separately reflect the landscapes f loc across corresponding interfaces, which are w ss for the poresubstance interface, w sl for the solid-liquid interface inside the substance, and w gr for the grain boundaries inside the solid, i.e.: In this regard, w ss can be formulated in a simple double-well form ð10Þ As mentioned, CðTÞ, DðTÞ and HðTÞ are the temperature-dependent barrier height which contains the parameters derived from the internal energy and configuration entropy: HðTÞ ¼ H pt À H cf ðT À T M Þ: Finally, the gradient term f grad is: f grad ðT; rq; fr/ i g; frg j gÞ ¼ Ts grad ðrq; fr/ i g; frg j gÞ ð11Þ Figure 2 presents the overall landscape of local free energy density f loc . Since f loc reflects the local thermodynamic stability of a certain phase at a specific temperature, it is anticipated to have multiple minima, including ( The heat term f ht tilts the multi-well due to the local temperature variation, manifesting the thermodynamic stability of a variable due to the change of local thermal conditions. To get a close view, we further separate f ht into two parts: the contribution from the heat absorbing/releasing due to local variation of temperature , and the contribution from These two contributions are also indicated in Fig. 2a and b. When local temperature meets the thermal conditions of melting (i.e., overheating, when T > T M ) and resolidification (i.e., undercooling, when T < T M ), vast tilting of the multi-well on the side of the liquid phase can be observed due to the additional contribution from f L . This indicates the thermodynamic driving force of the phase transition. On the other hand, for every temperature variation, grains with the same temperature should still maintain equality in local thermodynamic stability (illustrated as the equal multi-well among fg j g) for ideal grain growth, until the temperature gradient further breaks the local equilibrium among the grains. 39 This is achieved by the interpolating function U ht . Constraints of Eqs. 1 and 2 are also presented as the lower bounds in Fig. 2 to maintain the shape of the f loc at the corresponding threephase junctions. Note that the substance constraint in Eq. 1 is satisfied when / l ¼ q À / s . Therefore, we can safely use only / to represent the solid substance and ðq À /Þ to represent liquid substance. The other constraint in Eq. 2 is fulfilled using the penalty method. In this regard, the constraint term pRĝ with a penalty coefficientp and the penalty sum Rĝ is added in the free energy density. Eventually, the explicit formulation of F in Eq. 5 becomes: þf grad ðT; rq; r/; frg j gÞ þpRĝdX: f loc ðT; q; /; fg j gÞ ¼ w ss ðT; qÞ þ w sl ðT; q; /Þ þ w gr ðT; /; fg j gÞ; ð12Þ with w ss ðT; qÞ ¼ CðTÞ q 2 ð1 À qÞ 2 h i ; w gr ðT; /; fg j gÞ ¼ DðTÞ Interpolating functions U ss , U s and U l can be simply formulated as: and for U ht , we simply adopt the formulation: To include the constraint from Eq. 2, the penalty sum is formulated as: with the exterior penalty functionsĝ / andĝ g j as: Model parameters A; B; G; C pt ; C cf ; D pt ; D cf ; H pt ; H cf , as well as the gradient constants (j q , j / and j g ), are obtained from given diffusive interface width and the experimentally measured interface energies (i.e., c exp gs , c exp gl , c exp sl and c exp gb for the pore-solid, poreliquid, solid-liquid interfaces, and grain boundary). Coefficient n is employed to favor the determination of model parameters by fitting the experimental results. 33 A coupled kinetic system for OPs, a melting flow dynamic as well as heat transfer is then derived under the non-isothermal thermodynamic framework shown above. The derivations are explicitly presented in the Supplementary Note 1. As a result, we obtain the nonlinear kinetic system for the model as: .
@e @T DT Dt þ @e @q Equations 16a and 16b govern the flow of the melts with the flow velocity vector u and hydrostatic stress p as the coupled variables. m is the dynamic viscosity, b is the body force, and . is the density. The Cauchy stress tensor in which the Korteweg stress tensor r K reflects the effects of the interface tensions on the fluids, i.e., the capillary effect. 26,43 Since in this case r K is also a function of T, r Á r K thereby indicates the possible driving force due to the inhomogeneity of interface tensions which is induced by on-site temperature gradients. This is also known as the thermocapillary effect (or Gibbs-Marangoni effect). 26,44,45 Apart from this, Eq. 16c presents the thermal evolution equation coupled with the transient microstructure (i.e., transient terms of OPs) and melt flow dynamics. Equations 16e and 16f result in the form of the Allen-Cahn equations. 46 Equation 16d results in the form of the Cahn-Hilliard equation 47 with a direct temperature gradient-induced diffusion term, which is also known as the thermophoresis effect (or Ludwig-Soret effect). 48 The derived kinetic system in Eqs. 16a-16f can be regarded as the combination of the NSCH and NSAC systems coupled with a transient temperature field and temperature gradient-related terms, such as thermocapillary and thermophoresis. In the following simulation, however, the thermocapillary and thermophoresis terms will be tentatively dropped due to the consideration of the computational stability and cost. Those effects will be explicitly covered in upcoming works.

Implementation of Finite Element Method
The model is numerically implemented via the finite element method (FEM) within the program ''NIsoS'' developed by authors based on MOOSE framework. 49 4-node quadrilateral Lagrangian elements are chosen to mesh the geometry. The Cahn-Hilliard equation is solved in a split way. 50,51 A transient solver with a preconditioned Jacobian-Free Newton-Krylov method and backward Euler algorithm has been employed to solve the nonisothermal phase-field problems. Adaptive meshing and time stepping schemes are used to reduce the computation costs. In addition, to stabilize the calculation of NSCH and NSAC systems, streamline-upwind Petrov-Gelerkin and pressure-stabilized Petrov-Galerkin methods are introduced associated with the weak forms of the Navier-Stokes equations. 52

Simulation Domain
In this work, we consider a 2D simulation domain of the powder bed mid-section, as shown in Fig. 3. Instead of performing 3D simulations directly, 2D mid-section scenarios can already provide sufficient microstructure information of a powder bed, i.e., particles with multiple sizes and various pores due to the particle packing. Although details of other sections are lost, it can still recapitulate multiple features observed during PBF processes, including melting/partial melting, re-solidification and sintering of particles/grains, melt flow as well as pores/ bubbles behaviors, the thermal profile and its evolution, and resultant grain structures, without having a vast computational cost and severe issues of numerical stabilization. This also helps us to easily verify underlying physical processes and their interactions. The 2D simulation domain has a thickness of H pb and a length of l pb . Particles inside the domain are generated with the random close packing procedure. Due to the uncertainty of the initial grain structure of a single particle, we simply treat each particle as a monocrystal with an unique random orientation following the reported simulation works. 13,14,34,39 With the help of a minimum coloring algorithm and a grain-tracking algorithm, 33,53 four g i are sufficient to uniquely represent all the particles/grains for these simulations. Two type of boundaries (i.e., the atmosphere and substrate boundaries) are considered as the combination of three types of boundary conditions (BCs). The mass closed BC allows no mass flux or flow on the boundary ð17Þ wheren is the normal vector of the boundary C, 0 is the null vector. Next, the heat convective BC allowed heat dissipation as heat convection with the atmosphere where h is the convective coefficient and T E is the environment temperature. Finally, the heat conductive BC allowed heat dissipation in a way of heat conduction through the substrate where H sub and k sub are the thickness and the thermal conductivity constant of the homogeneous substrate, respectively. The bottom of the substrate is fixed at a temperature T P , the pre-heating temperature.
To distinguish properties on substance (ss), atmosphere (at), surface (sf) and grain boundary (gb), interpolating functions U ss , U at , U sf and U gb are utilized, which obtain unity only in the corresponding region, to construct overall properties from individual ones. They can be simply formulated as: ð20Þ Assuming the isotropic diffusion and grain boundary migration, the effective value of mobility L g and M path through possible paths (path ¼ ss, at, sf and gb) are adopted correspondingly from the selfdiffusivities D eff path and grain boundary mobility G eff gb ; 39,54,55 i.e.: Then, the diffusive mobility M on the powder bed can be formulated as 33,39 : Similarly, for isotropic heat conductivity k and density . we distinguish them on the substance (ss) and atmosphere/pores (at): and for viscosity m, the solid (s), melt (l) and atmosphere/pores (at): ð24Þ As mentioned, the thermal effect is equivalently treated as an internal heat source term q v moving with a velocity v: in which P 0 is the nominal beam power reaching the surface of the powder bed, b is the attenuation coefficient. (x, z) is an arbitrary point inside the simulation domain, while ðx v ; z v Þ is the moving center of the beam following the trajectory of the surface. Equation 25 can be regarded as the midsection of a volumetric heat source with a 2D Gaussian distribution along the surface and a decay along the depth 56,57 (also illustrated in Fig. 3), which still follows the scenario of this work. R bm and f are the characteristic radius and penetration depth of the beam. Note that parameter P is utilized to adjust the concentration of the deposited power inside the circular beam spot, e.g., when P ¼ 4:6 there is 99% of the P 0 concentrated inside the circular beam spot characteristic by R bm . f=2 indicates the half-decay depth of the deposited energy.

Parameter Normalization
Following the conventions of the NSCH and NSAC systems, dimensionless quantities, including the Péclet number (Pe), Reynolds number (Re), Cahn number (Ch), Weber number (We) and Froude number (Fr) are defined to characterize the scale of diffusion, viscosity, interface, capillary and inertial force comparing to the scale of flow, respectively: where v and . are the characteristic velocity and density of the flow, while c is the characteristic surface tension. m is dynamic viscosity and M is the isotropic diffusivity.
l is the characteristic diameter of the viscous flow, while ' is the characteristic width of the diffusive interface, and b is the inertial force (principally the gravity). For convenience, we use a set of simpler reference quantities to re-define those characteristic quantities, including the reference (melting) temperature T M , the reference length scale l, and the time scale t. The reference energy density C T M pt ¼ j T M q T M = l 2 which is the model parameter obtained at the reference temperature T M . j T M q is the gradient model parameter at a reference temperature T M . Then, we relate all of the characteristic quantities to the references as: ð27Þ Dimensionless quantities in Eq. 26 can be thereby modified as: note here that we multiply the Cahn number by the Weber number to obtain 1. Also, since multiple diffusion paths are considered, the Péclet number of each diffusion path Pe diff path should be specifically calculated. Similarly, as for the NSAC system, we define another two Péclet-type numbers Pe sl and Pe gb to characterize the mobility of the solid-liquid phase transition as well as grain growth, respectively Dimensionless forms of other quantities are given in Table I. Spatial and time derivatives are also normalized with respect to t and l, respectively. Dropping the tildes for simplicity, we eventually obtain the dimensionless forms of the kinetic equations in Eqs. 16a-16f as:

RESULTS AND DISCUSSION
The simulation results of the single scan of stainless steel 316L (SS316L) in an argon atmosphere are presented here. SS316L belongs to the family of austenite stainless steels. For this type of stainless steel, no solid structural transition is anticipated. 66 The reference length scale is set as l ¼ 1 lm, and the time scale as t ¼ 1 ls. Characteristic radius and depth of the beam are set as R bm ¼ 10 lm and f ¼ 80 lm. The attenuation coefficients are set as b ¼ 0:50. Scan speed v is simply defined as v ¼ jvj. The melting point of SS316L is set as T M ¼ 1700 K. The temperature of the powder bed is initialized as Tj t¼0 ¼ 0:4T M ¼ 680 K , which is also the value of the pre-heating temperature T P and the environment temperature T E . Utilized material properties are shown in Table II and the diffusive interface width of the grain boundary reads as . Due to the difficulties in obtaining the experimental value as well as the temperature-dependent trend of c exp gl and c exp sl , we assume the value of H pt , H cf and j / as shown in Table III. We set l diff ¼ 2 lm, parameters A; B; C pt ; C cf ; D pt ; D cf ; j g and j q are fitted from the experimental surface and grain boundary energy c exp gs and c exp gb through the method we proposed in our latest work. 33 Figure 4 shows the chronological progress of the single scan and highlights the multiple phenomena. A powder bed with H pb ¼ 100 lm and l pb ¼ 500 lm is utilized. The thickness of the substrate is H sub ¼ 200 lm. Due to the beam scan, the local temperature of the powder bed rapidly increases. A higher temperature gradient can be observed at the moving front of the beam than at the tail. Once the

Model quantities
local temperature exceeds the melting point T M (i.e., the region is overheated), the melt pool is formed. Inside the melt pool, different kinds of melt flow can be observed. Due to the existence of the gravitational Froude term, melts at the convex position tend to flow down, causing a flat flow around the surface (Fig. 4c and d) and floating pores (inset in Fig. 4b). Due to the Allen-Cahn-governed   phase transition (with the mobility L / ), we note that the melting process does not occur instantly inside the overheating region. As a result, partially melted particles/grains can be observed (Fig. 4a-d). Correspondingly, undercooled melts also exist (the melt phase in the region where T < T M , i.e., the undercooling region) where the re-solidification of grains can be observed. It is anticipated that, when the resolidification rate is approximately equal to the scan speed, a dragged tail of melt pool forms, behind which the columnar grains remains. Pores/bubbles which have not yet emerged to the surface are thereby trapped. Outside the melt pool, the local temperature is still high enough to activate diffusion and thus induces necking between adjacent particles ( Fig. 4a and b), presenting the feature of sintering. Apart from the captured phenomena, Fig. 4 also presents the interactions between the local morphology, temperature field, and kinetics such as the melting, re-solidification and grain growth. It is obvious that the temperature field (i.e. isotherms) follows the morphology of the powder bed, showing the locally concentrated temperature gradient around the necks among the powders. Correspondingly, those kinetics interact mutually with the local morphology and the thermal profile. For instance, the melt flows flatten the surface and enlarge the melt pool. Since the melt flow transports not only mass but also energy, it thus influences the local thermal profile. Similarly, re-solidification of the grains does not just follow the local temperature field, but it also creates a grain structure which is distinct from the one before the beam scans, resulting in a different local environment for heat transfer at the tail of the melt pool compared to the front. We further perform a series of simulations with varying beam power P 0 and scan speed v. In-process microstructures with melt pools are shown in Figs. 5 and 6 for increasing beam power (300 W-600 W, 2000 mm s À1 ) and decreasing scan speed (400 W, 3000 mm s À1 -1250 mm s À1 ) , respectively. The powder bed has the dimensions of H pb ¼ 100 lm and l pb ¼ 1000 lm. The thickness of the substrate is H sub ¼ 200 lm. The beam scans start at the horizontal position x ¼ 200 lm on the surface of the powder bed and end at x ¼ 900 lm. The melt pools generated with varying geometries corresponding to the different P 0 and v, move along with the beam scan. The dragged tails form for all cases, indicating the re-solidification has reached an equilibrium with the scan-induced melting (i.e., the resolidification rate at the rear of melt pool equals the scan speed). A few interesting observations can be noted. Firstly, the re-solidification front follows the local profile of the temperature field. More exactly, grains are re-solidified along the local cooling directions, i.e., the negative gradient of the temperature field. Secondly, the melt pool size varies differently from the increasing power than from the decreasing scan speed. Figure 5 demonstrates that when P 0 increases, both the length and the depth of the melt pool increase. On the the other hand, with decreasing scan speed, only the depth of the melt pool increases monotonically, while the length of the melt pool first increases and then decreases when v becomes small enough, as shown in Fig. 6. At the same time, the temperature gradient at the melting front becomes much higher. If we consider the linear power density as P 0 =v, it demonstrates the possible length reduction of the melt pool for a case with high linear power density. A possible reason is the enhancement of heat dissipation induced by the on-site morphological changes. Note that a less porous morphology is formed, which is preferred for heat conduction according to the implementation of k in Eq. 23. Therefore, in the cases with a low scan speed, although the linear power density is high enough to generate a melting pool with considerable depth, the rate of heat dissipation is competitive to the beam scan, resulting in the reduction in length of the melt pool. Apart from these, undercooled melts of all cases, regarded as the region from the isotherm of T ¼ T M to the dragged tail, are observed to have the same undercooling degree. This is due to the implementation of a constant mobility of phase transition L / (or the Péclet-type number Pe sl ). In a near-realistic system, the relationship between L / (Pe sl ) and undercooling degree should be properly formulated.
Melt pool geometries, varying with beam power and scan speed, significantly influence the morphology as well as the porosity/densification of the resultant microstructure, especially the morphology of the pores and grains. To systematically discuss such influences, we present the densification and corresponding microstructure of the PBF-processed powder bed as a distribution with different beam powers and scan speeds in Fig. 7. Densification factoris defined as: where e is the porosity of the final microstructure, while e 0 and e min represent the initial and minimum porosity achieved through the processes, respectively. In this work, e min ¼ 0. According to the morphologies, This map can be roughly divided into four regions, as shown in Fig. 7a. Cases located in the wide intermediate region (indicated by -¼ 10% $ 95%) result in the coexistence of columnar grains and irregular-shaped pores, which have been observed in the experimental works. 67,68 The formation of the irregular-shaped pores can be explained by the competition between the bubble dynamics inside a melt pool and the solidification. We track two sets of the pores on the processed powder bed, labeled as S 1 and S 2 . As shown in Figs. 5b-e and 6b-e, pores behave as floating bubbles, driven by the hydrostatic pressure due to the existence of the Froude term. During this process, large bubbles (S 2 ) deform, while small ones (S 1 ) reshape themselves into rounds or ellipses. When re-solidification occurs in the melt pools with relatively small sizes (Figs. 5b, c and 6c, d), large bubbles are trapped as irregular-shaped pores, yet their shapes are more or less similar to the original tunnel-like ones (S 2 in Fig. 7b 3 -b 6 ). Small bubbles may also be trapped as is, though some of them around the surface find it difficult to emerge. As a result, there are concentrated small pores around the surface (S 1 in Fig. 7b 3 -b 4 ). When the melt pools are relatively large (Figs. 5e and 6a), large pores may be significantly split into small ones and quickly ejected from the pools. Once trapped, they also present highly deformed/split shapes (S 2 in Fig. 7b 1 -b 2 ). In addition, the columnar grains are also observed to have different heights corresponding to the melt pool size that has actually been achieved. When increasing the scan speed while Simulation results on PBF processing of the SS316L powder bed with a constant v ¼ 2000 mm s À1 with varying P 0 , when the beam is at x ¼ 700 lm on the surface of powder bed. Tracked pore sets S 1 (dash) and S 2 (dash-dot) are also marked. Fig. 6. Simulation results on PBF processing of the SS316L powder bed with a constant P 0 ¼ 400 W with varying v, when the beam is at x ¼ 700 lm on the surface of powder bed. Tracked pore sets S 1 (dash) and S 2 (dash-dot) are also marked.
fixing the beam power (compare Fig. 7b 1 , b 3 and b 6 ), the depth of the melt pool reduces, resulting in grains with less height. This works in a similar fasion when decreasing the beam power while fixing the scan speed (compare Fig. 7b 2 , b 3 and b 5 ). Cases of Fig. 7b 2 and b 3 locate at almost the same isoline of -, they present similar morphology of columnar grains and pores. Apart from this, packed tilted columnars with almost no irregular-shaped pores (Fig. 7c 1 -c 3 ) can be observed in cases located on the higher-left region (indicated by -! 95%Þ), which also presents the resemblance to the experimental observations on the SS316L. 69 In those cases, significant melt pools with depths larger than the thickness of the powder bed are generated, and almost all pores are split and ejected from the melt pool. The number of residual pores is limited, and their shape is mostly round. However, they are found in the grain interior, at the grain boundary and at grain triple junctions. Resolidified grains follow the local cooling directions and eventually form a tilted column shape. Note that grains at the end of scans tilt in the other direction from the beginning ones due to the reversal of the local cooling direction. It is worth remarking that intense vaporization can be expected in this region when the local temperature is over the boiling point (possibly caused by either sufficiently high beam power or sufficiently low scan speed). In this sense, the recoil pressure induced by vaporization applies additional force on the liquid-gas surface and creates a depression of the melt pool. As a result, there is a possible drop of the densification due to trapped extra pores after refilling of the depression. 8 Since the related effects of the vaporization are not yet contained in this work, details about this densification drop are thereby lost.
As for cases located in the narrow strip-shaped region (indicated by -¼ 4% $ 10%), the melt pools are rather limited and even fragmented, as shown in Fig. 7d 1 . Regarding melt pool fragmentation, there are different explanations, including flow instability or powder morphology. 5,10,70 Since effects such as wetting and Marangoni flows have not yet been implemented in this work, the simulated melt pool fragmentation here should be attributed to the local stochastics of the powder bed. Enhanced heat dissipation induced by the on-site morphological changes is also a possible reason, causing the inhomogeneous heat loss from the melt pool. The resultant microstructure (Fig. 7d 2 ) inside this region presents similarities to those observed from liquid-state sintering where the particles are sintered with the existence of the liquid phase. 33 Finally, cases located in the right-bottom region (indicated by -4%) do not have melting phenomena, and barely present features observed from solid-state sintering except for necking.

CONCLUSION
In this work, we propose a unified non-isothermal phase-field model for the heat-melt-microstructure-coupled processes during PBF, which consists of a non-isothermal free energy density functional considering effects from pore/liquid/solid phases as well as the grains in the solid phase, and a framework deriving simultaneously the interactive kinetics for thermal evolution, melt flow dynamics and microstructure evolution. Using the proposed model, finite element simulations on PBF of SS316L were carried out. The results clearly demonstrate the promising capability of the model. The simulation results are briefly summarized in the following: 1. The proposed model is able to recapitulate phenomena during PBF, such as high-gradient temperature field, melt pool with various geometry, melt flow as well as pores/bubbles behaviors, partially melted particles, re-solidification and sintering of particles/grains, and resultant columnar grains with irregular-shaped pores. It also presents the interactions between transient temperature field, melt flow dynamics and grain structures, e.g., the enhanced heat dissipation due to the formation of on-site morphology with less porosity, and tilted re-solidified grains following the local cooling directions. 2. Simulation results on the melt pool indicate strong interaction between the melt pool and the powder bed morphology. Parameter studies imply that the melt pool length varies in different manners with respect to the increasing power and to the decreasing scan speed. The length increases monotonically with increasing power, but for decreasing speed, it first increases and then starts to decrease after the speed becomes small enough. This may be attributed to the competition between heat dissipation, enhanced by on-site morphological changes, and scan-induced melting. 3. A densification map with respect to the beam power and the scan speed was obtained. The isolines of the densification factor are not following the isolines of the specific line energy, indicating again that the beam power and the scan speed should be addressed as independent process parameters. Four regions of the parameter combinations are identified based on the resultant microstructure and pore features. Inprocess simulation results allow us to trace back the formation mechanism of the various types of irregular-shaped pores.
The current version of the proposed model should be further extended in the near future for different aspects, e.g., implementation of capillary/thermocapillary effects (Korteweg stress term), and the thermophoresis effect (direct temperature gradientdriven mass transfer). In order to consider the effects from the recoil pressure induced by vaporization, the pressure balance should be explicitly formulated at the liquid-gas surface via asymptotic analysis. In addition, quantitative description of the model should be further derived and examined. from the HHLR, Technische Universitä t Darmstadt.

ELECTRONIC SUPPLEMENTARY MATERIAL
The online version of this article (https://doi.org/ 10.1007/s11837-019-03982-y) contains supplementary material, which is available to authorized users.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons li-cence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecom mons.org/licenses/by/4.0/.