Three-dimensional non-isothermal phase-field modeling of microstructure evolution during selective laser sintering

Predicting the microstructure during selective laser sintering (SLS) is of great interests, which can compliment the current time and cost expensive trial-and-error principle with an efficient computational design tool. However, it still remains a great challenge to simulate the microstructure evolution during SLS due to the complex underlying physical phenomena. In this work, we present a three-dimensional finite element phase-field simulation of the SLS single scan, and revealed the process-microstructure relation. We use a thermodynamically consistent non-isothermal phase-field model including various physics (i.e. partial melting, pore structure evolution, diffusion, grain boundary migration, and coupled heat transfer), and interaction of powder bed and laser power absorption. The initial powder bed is generated by the discrete element method. Moreover, we present in the manuscript a novel algorithm analogy to minimum coloring problem and managed to simulate a system of 200 grains with grain tracking using as low as 8 non-conserved order parameters. The developed model is shown to capture interesting phenomena which are not accessible to the conventional isothermal model. Specifically, applying the model to SLS of the stainless steel 316L powder, we identify the influences of laser power and scanning speed on microstructural indicators, including the porosity, surface morphology, temperature profile, grain geometry, and densification. We further validate the first-order kinetics during the porosity evolution, and demonstrate the applicability of the developed model in predicting the linkage of densification factor to the specific energy input during SLS.


Introduction
Selective laser sintering (SLS) is a typical additive manufacturing process meant for rapid prototyping and tooling [1][2][3][4][5]. During SLS, a desired geometry is built by sequentially layer-by-layer powder spreading and subsequent laser scanning, whereby the photonic energy is transformed into heat by absorption [6,7]. To be distinguished with the other laser-based powder bed additive manufacturing known as selective laser melting (SLM), there is no significant melting phenomenon during SLS. Whereas the temperature is sufficiently high for particles to bind together through sorts of mechanisms, leading to the products with relatively high porosity [5,8,9]. Because of these characteristics, SLS has been applied for the industrial production of individually designed components made of organic polymers, ceramics, and metallic alloys which have relatively high melting or transition temperature [1,5,10,11]. It also shows the possibility to produce the porous biomaterials, especially medical scaffold and bones [12][13][14][15], and functional materials [10,16,17].
Although the procedure of SLS is conceptually simple, the underlying physics is complex and covers a broad range of time and length scales. Since its birth [18], identification and comprehension of those phenomena during SLS and their effects on the final product are crucial for the successful manufacturing, which highly relies on trial-and-error and variation across the grain boundary is negligible). When ρ = 0, no grain is present. This profile of OPs leads to the constraint (1 − ρ) + α η α = 1. Their temporal evolution indicates the changes of the surface and grain boundaries, and thus represents the microstructure evolution during the SLS process.
The non-isothermal formulations following our latest work in Ref. [39] are extended with the contribution from partial melting. The total free energy density f (T, ρ, {η α }), derived through the Legendre transform of the internal energy and entropy [41], is formulated as where f ht (T ) = T e ht (T )d 1 T , The eight model parameters (A, B, C pt , C cf , D pt , D cf , κ η and κ ρ ) in Eq. (1) are obtained from the experimentally measured temperature-dependent surface and grain-boundary energies (γ sf and γ gb ) and the width of the grain boundary. Fitting coefficient ξ is utilized to help fitting with the experimental result (see Supplementary Information). e ht is the gain (or loss) of the internal energy density due to the heat transfer, which can be formulated as (note that we have set the internal energy of the atmosphere/pores at melting point T M of the system as zero) In Eq. (2), c is the volumetric specific heat. The subscript 'bk' and 'at' denote the substance materials and atmosphere/pores, respectively. The latent heat during the solid-liquid phase transition is simply represented by L Φ(τ ), where the sufficiently smooth interpolation function Φ(τ ) takes one when τ = T /T M → 1 and to zero when τ → 0. The temporal evolution of ρ is derived from the mass conservation, i.e.
where j diff and j melt represent the mass flux contribution from the diffusion and partial melting, respectively. Here we simply assume that the partial melting is locally restricted to individual particles, and the melt flow is driven by the capillary pressure. The driven force of the melts is described by change of surfacial curvature A sf and energy γ sf , i.e. j melt ∼ −Φ(τ )M melt · ∇(A sf γ sf ). According to Young-Laplace equation in which the on-site chemical potential is proportional to A sf γ sf , we can thereby formulate j melt into a diffusion-like form. j diff can be obtained directly from Fick's law. In this way, Eq. (3) can be rearranged into the form of the Cahn-Hilliard equatioṅ Here the mobility tensor M is formulated in the fashion to consider contributions not only from mass transfer paths through bulk (bk), atmosphere (at), surface (sf) and grain boundary (gb) [34,40,42,43], but also from the gain due to the partial melting, i.e.
In this regard the contribution of the partial melting is treated more like an enhanced surface diffusion when τ → 1.
To represent such contribution more properly, a j melt obtained by coupling with fluid dynamic equations (e.g. the Navier-Stokes Equations) should be worked out in the future. On the other hand, evolution of {η α } is governed by Allen-Cahn equation with the corresponding mobility tensor The laser scanning process is modeled by moving the Gaussian distributed heat source on the powder-bed surface at a speed of v. The kinematic equation of temperature is formulated as where the thermal conductivity k adopts spatial distribution as M in Eq. (5) i.e.
q(r) represents the volumetric energy deposition due to the radiative energy flux of the laser and can be formulated as [6,7] q(r) = −βP 0 p xy (x, y)p z (α, λ, z), (9) in which P 0 is the nominal laser power reaching the surface of powder bed, p xy (x, y) is the 2D Gaussian distribution. p z (α, λ, z) is a penetration function, which takes the form proposed by Gusarov et al. [6]. The hemispherical reflectivity α indicates the influence of the powder material. The attenuation coefficient β and the optical thickness λ reflect the influence from the powder bed structure. For a loosely packed powder bed with particles of distributed diameters d, a thickness of H pb , and a porosity of ε, the effective value of these two parameters are calculated as Since the details inside the substrate is trivial in this work. we set a heat conduction boundary condition (BC) to replace the mesh of the a substrate (Fig. 1a), which has a similar fashion to the third-type BC of heat transfer problems, i.e.
where n is the normal vector of the boundary Γ sub . H sb and H sb are the thickness and the heat conductivity constant of the homogeneous substrate, respectively. The bottom of the substrate is set with a fixed temperature T P as the pre-heating temperature.

Optimized profile of order parameters
It is anticipated that a general powder bed would have a spatial and temporal conserved OP ρ(r, t) in combination with the multiple non-conserved OPs {η α } = {η 1 (r, t), η 2 (r, t), . . . η N (r, t)}. The number of non-conserved OPs N is required to be the same as the number of the particles/grains to statistically retain the uniqueness of each crystalline orientation [44,45]. However, such computation is expensive due to the large amount of variables, and inefficient since only a few of variables have nonzero value at any point in the domain. To solve those problems, methods that can drastically reduce the computational cost while retaining the uniqueness of the particles/grains have been proposed in recent decade [46][47][48]. The fundamental idea is to assign the same OP to grains which are sufficiently spaced, and remap the OPs when certain grains with an identical OP tend to have coalesced. This idea makes it possible to use less amount of non-conserved OPs to simulate the evolution of the polycrystalline structure. Taking the example of the grain tracking algorithm proposed by Permann et al. in [48], it is able to simulate more than 1,000 grains, and only a minimum of 8 non-conserved OPs in 2D and 28 in 3D are required for the grain growth problems. Due to features of SLS, like on-site heating and rapid cooling, global and long-term grain coalescence is absent, in contrast to the grain growth process. This means usually no severe changes of the profile of non-conserved OPs take place, making it possible to further reduce the computational cost by optimizing the profile of OPs.
In this work, we translate the problem of OPs' assignment into a classical minimum coloring problem (MCP). By solving the MCP with the Welsh-Powell algorithm [49], we obtain the optimized profile with a minimum number of non-conserved OPs. The work flow of the process is presented in Fig. S1a (Supplementary Information). In the beginning, an adjacency matrix of the powder bed (Fig. 1c), demonstrating the "adjacent neighbors" of each powder, is generated according to the criterion of adjacency. If two particles/grains which are spaced less than this critical distance, they can be considered as "adjacent neighbors". Here the distance between the largest powder and its nextnearest neighbor in hexagonal close packing (inset of in Fig. 1c) is considered as the criterion of adjacency. Then, the map G(V, E) is established based on the adjacency matrix, which consists the unique particle ν i (r, R) as the vertices V{ν i } and corresponding degree of adjacency E{ i }, i.e. there are i particles are adjacent to the one indexed as ν i (Fig. 1d). Once the map is established, the algorithm will iteratively assign the available non-conserved OPs to non-adjacent and non-assigned particles till all of them have been assigned with OPs. The particles indices remain identical during this iterative assignment, and can be later translated by the grain tracking algorithm ( Fig. 1e and  1f). In this work, for the powder bed with particle diameters ranging from 20 to 50 µm, 8 non-conserved OPs are sufficient to represent about 200 particles/grains. For another simulation attempt where a loosely packed powder bed with about 400 uniformly sized particles is utilized, only 6 non-conserved OPs are sufficient. In this way, the number of non-conserved OPs is remarkably reduced, and the computation resource is efficiently saved.

Results
Here we present the simulation results on the microstructure and surface morphology evolution by using the experimental data of the type 316L stainless steel (SS316L) in the inert atmosphere. Assuming the isotropic diffusion and grain boundary migration, the mobilities L and M are formulated as scalars and adopted from the self-diffusivity D eff path through possible paths (path = bk, vp, sf and gb) and grain boundary mobility G eff gb [36,39,50], i.e.    Fig. 2a presents the microstructure evolution obtained at a laser power of 20 W and a scanning speed of 100 mm/s. For realizing the SLS process, a laser beam then scans the powder layer straightly and binds the free particles together simultaneously to form a continuous but porous morphology. The SLS induced physical phenomena around the laser beam spot are shown in Fig. 2b. Due to the laser scanning, temperature of the powder bed obviously increases, as shown in Fig. 2a. Some regions even have the temperature over the melting point T M , creating the thermodynamic condition for the physical processes such as diffusion and partial melting. The temperature is inhomogeneously distributed and only the particles around the laser beam spot are partially melted. In the region with T > T M , the reduction of surface energy (or surface capillary pressure) makes the melts flow from the convex to the concave region. This kind of melt flow results in a continuous region with coarsened grains once cooled down to T < T M . In the region with a maximum temperature lower than T M , no melting occurs. But the temperature is still high enough to activate diffusion and thus induce necking among adjacent particles. In addition, from the isotherms in Fig. 2a, we can estimate the local temperature gradient around the partial melting region and the pore region as high as 50 K/µm and 100 K/µm, respectively. Such large-gradient temperature field is essential to the non-isothermal behaviors of the SLS, which could influence the grain coalescence and eventually affect the microstructure evolution [39].
Laser power and scanning speed are important SLS parameters through which one can effectively control the properties (e.g. porosity and grain size) of the processed components. Under a constant laser power of 5 W, we present the simulated surface morphologies along with the corresponding experimental observations in Fig. 2c and 2d for a scanning speed of 20 and 5 mm/s, respectively. It is obvious that changing the scanning speed has significant effects on the attained morphology. For the scanning speed of 20 mm/s in Fig. 2c, powder bed shows a poor binding and the necking among adjacent particles is very weak. For the scanning speed of 5 mm/s in Fig. 2d, on the other hand, a continuous piece form after the SLS processing, indicating the better binding of particles.
Furthermore, we investigate the influence of scanning speed and laser power by performing a series of simulations. As shown in Fig. 3, the morphology and porosity of the SLS processed powder bed are presented with different laser power and scanning speed. The porosity map can be roughly divided into two regions by the dot-dash line as shown in Fig. 3a. The lower-right region shows less bound particles, while the upper-left region shows more continuous region (or even fully melted as shown in Fig.3b 4 ). When decreasing the scanning speed while fixing the laser power, more particles are bound to create more continuous pieces, i.e. less porosity will be achieved in the final components (comparing Fig. 3b 4 -3b 7 ; 3b 8 -3b 10 ; and 3b 3 , 3b 12 -3b 11 ). It is similar for the case when scanning speed is fixed and the laser power is increased (comparing Fig. 3b 1 -3b 4 ; 3b 5 , 3b 8 , 3b 12 ; and 3b 6 , 3b 10 -3b 11 ).
In the next section, we would explicitly discuss how laser power and scanning speed influence other features of the powder bed during SLS, including the temperature profile, particles/grains geometry evolution, and microstructure densification.

Temperature profile
Temperature profile in SLS is directly related to the processing parameters. It is an important indicator for judging whether certain physical process (e.g. partial melting) could occur during SLS. However, due to the difficulties in modeling the powder bed topography and the underlying physics strongly coupled with the heat transfer, simplifications are often used in simulating the temperature profile during powder bed based AM. For instance, the powder bed is frequently assumed to be a homogeneous bulk material with effective thermal properties [6,24,27,51]. There is also study in the particle scale to simulate the thermal evolution in a 3D multi-layer powder stacking model [26], but the 'particles' are kept in the cubic shape and packed regularly. In our work here, the powder bed generated by DEM consists of randomly packed particles. Moreover, the substance field ρ is introduced to assign different thermal properties to the substance and pores separately. The microstructure descriptor ρ and {η α } are also fully coupled with temperature. In this way, we can not only model the temperature field in the particle scale in a way close to the realistic setup, but also predict the thermally coupled microstructure evolution. Fig. 4a shows the profile of volumetric energy deposition from the laser beam, in which the half-maximum (HM) region is indicated by the dash line (diameter of FWHM), and the nominal beam spot by the dash-dot line (diameter of 2 × FWHM). Fig. 4b presents the top view and the section view of the temperature profile at the point when the laser beam has passed 350 µm on the powder bed with different processing parameters. It generally shows an inhomogeneous distribution of isotherms on the powder bed. In the front of the moving beam spot, isotherms are very dense and become sparse once the beam spot has passed away, indicating rapid heating up and a relatively slow cooling down. In this regard, increasing the laser power and decreasing the scanning speed can both improve the heat accumulation around the beam spot. In the case of P 0 = 20 W and v = 100 mm/s, we can only see some particles are partially melted. When increasing the laser power (e.g. P 0 = 30 W and v = 100 mm/s in Fig. 4b) or decreasing the scanning speed (e.g. P 0 = 20 W and v = 80 mm/s in Fig. 4b), the melting becomes more noticeable and the partial melting is even extended to full melting. In these cases mechanisms of melting and solidification dominate. On the other hand, decreasing the laser power (e.g. P 0 = 15 W and v = 100 mm/s in Fig. 4b) or increasing the scanning speed (e.g. P 0 = 20 W and v = 150 mm/s in Fig. 4b) can effectively reduce the melting, and in these cases the mechanism of solid-state sintering is dominant.
To systematically analyze the influence of processing parameters on the temperature profile and classify the possible dominant mechanism involved in SLS processing, we map the width (normalized with FWHM) and the depth (normalized with H pb ) of the melt with processing parameters, as shown in Fig. 4d and 4e. To help the discussion, each map is divided into four regions and the processing parameters inside a certain region are corresponding to similar microstructure features. For the normalized melting map (Fig. 4d) we have: (w1) no melting; (w2) extremely localized melting within the HM region (diameter of FWHM); (h3) wider melting region, yet still localized within the beam spot (diameter of 2 × FWHM); (w4) large-area melting. For the normalized melting depth map (Fig. 4e) we have: (h1) no melting; (h2) partial melting of the particles; (h3) some particles are fully melted and the melt depth exceeds half of the powder bed thickness; (h4) the powder bed is fully melted. Therefore, there is only dominant solid-state sintering in the overlapped region of (w1) and (h1). The dominant melting and solidification occur in regions (w4) and (h4). In the regions (w2), (w3) and (h2) and (h3), particles coexist with melts and thus the liquid-state sintering should be dominant. However, in the overlapped region of (w2) and (h2), the partial melting is restricted to the individual particle and then quickly regrows into the grain after cooling down. This process is closer to solid-state sintering. In other regions, there are more melts mixed with the unmelted particles, leading to typical liquid-state sintering.

Particles/grains geometry evolution
Due to the complex temperature profile, particles/grains during SLS may undergo different physical processes depending on the local thermodynamic conditions. Some particles around the beam spot may suffer from the partial melting. Others may still be sintered together since the temperature is sufficiently high to activate the diffusion and the grain boundary migration. Sorts of physics eventually lead to the evolution in particles/grains geometry evolution, including the change of grain size and grain shape. Such evolution, however, differs hugely from that during the conventional sintering where the temperature is firmly controlled. Since the high temperature gradient exists (there is in maximum 50 K/µm within the particles and 100 K/µm within the pores around the sintered neck), the diffusivity (or grain boundary mobility) would have a very large difference between the hot and the cool ends of particles around the partial melting region. For instance, the surface diffusivity at the hottest end of a partial melted particle is about 100 times larger than on at the coolest end in this work. According to our latest investigation in [39], highly non-isothermal condition can also lead to the difference in thermodynamic stability of grains. With the help of the grain tracking algorithm, we can track the particle/grain located at any position in the powder bed and investigate the geometric evolution during the process. Fig. 5a present the geometry evolution of particles/grains of the same processing parameters of P 0 = 20 W and v = 100 mm/s. Here the volume and the average temperature of the i-th particle/grain are calculated as where δ(Ω i ) is the delta function which takes unity inside the domain Ω i of the i-th particle/grain. Ω represents the simulation domain. The volume variation ∆V is thereby calculated as the absolute difference between the current volume V and the initial volume V 0 . The volume variation is normalized with respect to their initial values. When the laser beam passes, the average temperatures rise to the peak, then gradually drop. The volume variation, however, starts to rise when the average temperature reaches 0.6T M . After reaching the maximum T avg , the rate of volume variation gradually drops, and the grain volume tend to be constant. Grains A, D and E (Fig. 5a) with the same y coordinate shows almost the same volume variation. At the end of the scanning (5000 µs), they are in a polyhedral shape with arched surfaces (g A ,g D and g E ), presenting the resemblance to the grains during the intermediate stage of the sintering. Grain A, B and C with the same z coordinate present different geometric. Among these three grains, grain B is located farthest from the melting region and shows smaller volume variation and rate. It eventually turns into a shape which is less polyhedral, but similar to the one formed by pure necking without subsequent packing (g B ).
Grain C is located closest to the center of the partial melting region, whose T avg firstly rises above T M then quickly drops. It eventually turns into a flat grain (g C ) with a smooth top surface (due to partial melting) and polyhedral bottom surfaces (due to packing with other grains) . Grain F is exactly located beneath the partial melting region. Its volume variation increases firstly and drops later. Its volume variation reaches the maximum, then drops after reaching maximum T avg . In the end, it also turns into a polyhedral shape (g F ).
In Fig. 5b we also present the volume variation and variation rate of the same particle/grain (grain A) with different processing parameters. It can be found that increasing the laser power and decreasing the scanning speed can increase the volume variation. For v = 100 mm/s, the grain is polyhedral at P 0 = 20 W (g 3 ), but is round with a top-surface necking (g 5 ) at P 0 = 10 W. For the cases of P 0 = 30 W, v = 100 mm/s and P 0 = 20 W, v = 80 mm/s, flat grains with a smooth top and a polyhedral bottom are obtained (g 1 and g 2 ). We can also observe the fluctuation of the volume variation in the case of P 0 = 30 W, v = 100 mm/s when T avg is approaching its maximum, which is close to the melting point. This may be due to the contribution from the partial melting term Φ(τ )M melt .

Microstructure densification
Particle/grain geometric evolution where the grains get coarser and more packed while the pores shrink and reshapes, leads to the microstructure densification, which can be characterized by the porosity evolution. In Ref. [19] Simchi proposed that the temporal evolution of the porosity during the direct laser sintering of metal follows the first order kinetics law asε = −kε, where k is defined as the sintering rate constant. In Fig. 6a we present the temporal evolution of ln(ε/ε 0 ) of the three selected segments and the whole powder bed with an initial porosity ε 0 under the processing parameters of P 0 = 20 W, v = 100 mm/s. The average temperature T avg of each segment is also calculated according to Eq. (13). These three segments are with the same geometry (length 100 µm, width 250 µm), but located at different positions in the powder bed. It can be seen that when the laser beam scans into a certain segment, T avg rises and ln(ε/ε 0 ) is almost linearly decreased there. Once the laser beam is about to leave the certain segment, the T avg reaches its maximum then start to drop, while ln(ε/ε 0 ) is slowly decreased. In contrast, ln(ε/ε 0 ) of the whole powder bed presents an approximately linear trend, with a fitted slope as k = 7 × 10 −5 s −1 , as shown in Fig. 6a. This demonstrates that Eq. (13) is valid in this model to describe the porosity evolution of the powder bed combining all contributions from each finite segments, which is also in agreement with the experimental measurement in Ref. [19]. According to Fig. 6b, we can see the final microstructure (t = 5000 µs) around the segment center shows the resemblance to the microstructure during the intermediate stage of sintering [52]. In this stage, grains have been already packed together while most of the open pores have been eliminated or rearranged into a roundly closed shape. Whereas the microstructure around the segment margin presents less packed but more necked grains surrounded by tunnel-like open pores, manifesting the characteristics during the initial stage of sintering [52].
In Fig. 6c we present the temporal evolution of ln(ε/ε 0 ) of the whole powder bed under different processing parameters. Segment II is selected to present the influence of processing parameters on temporal T avg , which shows a similar trend as in Fig. 6a. The approximately linear trend of ln(ε/ε 0 ) is still obvious. Both increasing the laser power and decreasing the scanning speed result in the increase of maximum T avg and the rate constant k. Final microstructures of the segment II under various processing parameters are shown in Fig. 6d. We set the microstructure of segment II of P 0 = 20 W, v = 100 mm/s as the reference and compare it with cases with other processing parameters.  Here we further discuss the relation between the processing parameters and the sintering rate as well as the final densification. We define the specific energy input Ψ and densification factor Ω as In Eq. (15), P and v is the laser power and the scanning speed of the laser beam, respectively. h and w is the thickness and width of the scan track, respectively. ε min is the minimum attainable porosity in a sintered part, which ranges between 0.02 and 0.3 depending on the material properties. To avoid the deviation due to the less densified margin of the powder bed, we choose a FWHM = 100 µm of the power density of the laser as the track width w, and P = 75.8%P 0 (the surface integral of the power density reaches 75.8% of P 0 ). h takes H pb = 65 µm and ε min = 3%. As shown in Fig. 7a, the relation k vs Ψ can be linearly fitted with a slope of 9.51 × 10 −3 . Similarly, linear fitting of − ln(1 − Ω) vs Ψ give the densification coefficient K = 18.97, as shown in Fig. 7b. It should be noted that K is related to the material properties and the particle diameter distribution in the powder bed [19]. The K value obtained here is in line with the experimental K = 12.6 in Refs. [53] and [19] where the SS316L particles are with diameter < 45 µm.
The agreement in K shows the applicability of the developed model in simulating the microstructure densification. We also present a detailed map of the densification factor Ω in Fig. 7c. Following the isoline of the specific energy input Ψ (dot-dash lines in Fig. 7c), increasing laser power/scan speed actually leads to the increase of the densification Ω. And when higher Ψ are kept, the faster increase of Ω can be found. It reveals another interesting fact that specific energy input may not uniquely identify the densification of the processed component during SLS. Similar pattern has been experimentally observed in Ref. [54].
To sum up, we have developed a 3D non-isothermal phase-field model for investigating the evolution of the microstructure during SLS. Through numerical simulations, we capture interesting phenomena during SLS, such as temperature field with high gradient, mass transfer through partial melting, diffusion and evaporation, particle/grain necking and coarsening which leads to the microstructure densification. We also reveal the influences of the processing parameters (i.e. laser power and scanning speed) on microstructural features, including the porosity, surface morphology, temperature profile, particle/grain geometry, and densification. We further present the feasibility of the first-order kinetics during the porosity evolution, and verify the linkage of the sintering rate constant and the densification factor to the specific energy input. It is expected to further investigate the influence from factors such as the scanning strategies, initial particle size distribution and shape during SLS processing in the near future.

Powder bed generation using discrete element method
To generate the initial powder bed, the discrete element method (DEM) simulator YADE is used [55]. The spheres are firstly created as the gaseous loose-packing cluster with no contacts between particles. The distribution of the particle diameter is shown in Fig. S1b (Supplementary Information). Then, a gravitational force is imposed on the particles to make them spread into a rectangular box. The iterations continue until the deviation from force equilibrium on the particles is below a certain threshold, i.e the particles are stationary. This process is shown in Fig. S1b (Supplementary Information). After the powder bed is generated, the center r and radius R of each unique particle is recorded as a vertex ν i (r, R), and added into the vertices list V{ν i } for further optimization.

Implementation of finite element method
The model is numerically implemented within the MOOSE framework by finite element method (FEM) [56]. 8node hexahedron Lagrangian elements are chosen to mesh the geometry. The Cahn-Hilliard equation in Eq. (4), which is a 4th order differential equation, is solved by splitting it into two 2nd order differential equations via the introduction of an additional coupling field µ [57][58][59]. Then, the weak forms of split Eq.(4), along with the Eqs. (6) and (7), are discretized following the Galerkin method. To solve those time-dependent partial differential equations, transient solver with preconditioned Jacobian-Free Newton-Krylov (PJFNK) method and backward Euler algorithm are employed. Adaptive meshing and time stepping schemes are used to reduce the computation costs. The constraint of the order parameters is fulfilled using penalty functions.

Parallel CPU computation
The large-scale parallel CPU computations for each simulation domain, which has degree-of-freedoms (DOFs) on the order of 10,000,000 for both nonlinear system and auxiliary system, are performed with 150 processors and 2 GByte RAM per processor based on OpenMPI. Each simulation consumes on the order 10,000 of core hours.

Data availability
The authors declare that the data supporting the findings of this study are available within the paper and its Supplementary Information files.

Acknowledgement
The     Table S1 gives the material properties of SS316L in the argon atmosphere. Table S2 gives the dimensionless forms of the quantities utilized in this work. The dimensionless forms are obtained by normalizing with respect to a set of reference quantities, including the reference (melting) temperature T M , the reference length scale λ, and the time scale τ . The reference energy density C TM pt = κ TM ρ T M /λ 2 which is the model parameter obtained at the reference temperature T M . κ TM ρ is the gradient model parameter at a reference temperature T M . Other quantities are given in Table S2. Spatial and time derivatives are also normalized with respect to τ and λ, respectively. α316L 0.7 [7] * Activation energy is obtained from [8] while the prefix factor is estimated as unity at T M after normalization. ** Estimated as 100M eff sf .

Supplementary Note 1: Material properties and normalization
Supplementary Note 2: Obtaining model parameters by fitting the temperature-dependent surface and interface energies In this section, we present the derivation of the surface and grain-boundary energies at equilibrium which is based on the previous works by [9], [10], and [11]. Here we firstly show the dependency of the model parameters from Eq. (S1) to (S7), then derive the explicit formulation of the surface and grain-boundary energies from Eq. (S8) to (S15). Finally, direct and indirect methods of determining eight model parameters are proposed in the rest of this appendix. Considering the profiles of order parameters across the surface of a grain, where ρ and only one η varies from a semi-finite atmosphere/pore phase (−∞) to a semi-finite grain phase (+∞), as shown in Fig. S4. As for the normal direction r of an arbitrary point on the surface, profiles of ρ and η should satisfy the boundary conditions Then the specific surface energy can be calculated by Cahn's approach [12,13] where At equilibrium, functional in Eq. (S2) maintains minimum at each temperature T , requiring ρ and η to satisfy the Euler-Lagrange equation In this model, the constraint of order parameters (1 − ρ) + α η α = 1 should hold in any region within the substance at any time. Assuming ρ and η adopt the same shape as shown in Fig. S4a, i.e. ρ(r) = η(r), from Eq. (S3) we can yield the following equation where Replacing every η by ρ, we obtain To make Eq. (S5) hold at any T and ρ, one can assume where Known that A + B = 1, strong constraint above can be also rewritten as Eq. (S6b) shows the relation between the model parameters A, B, C pt , C cf , D pt , D cf , κ ρ and κ η . We can furthermore yield the expression γ sf for the specific surface energy. After replacing ρ by η in Eq. (S2), one can get the following Euler-Lagrange equation.  By integrating both sides of Eq. (S13) one can get f (T, 1, {η α , η β }) − T κ η |∇ r η| 2 = C 2 , (S14) and the arbitrary constant C 2 equals zero considering the boundary conditions in Eq. (S10). Substituting Eq. (S14) into Eq. (S11), we eventually obtain the explicit formulation of the specific grain-boundary energy at equilibrium (S15) We regard Eq. (S9) and (S15) as the relations which link the model parameters A, B, C pt , C cf , D pt , D cf , κ η and κ ρ to the measurable, temperature-dependent material properties γ sf (T ) and γ gb (T ), which are usually linearly fitted in practical measurement. Due to the dependency in Eq. S6b, there are only four independent quantities in these eight model parameters, i.e. D pt , D cf , κ ρ , κ η (or C pt , C cf , κ ρ , κ η ). Here we present a method utilizing fitting parameters to fit the temperature-dependent γ sf and γ gb to avoid the direct solving by formulating 4 equations based on Eqs. (S9) and (S15). This method firstly relates the κ ρ , κ η and C pt , D pt to the experimental values of the surface and grain-boundary energies approaching the melting temperature T M , i.e. γ TM sf = 2 where C pt and D pt are related according to Eq. (S6b). In the phase-field model, the grain boundary is a diffusive interface, and its approximate width λ gb at T 0 can be expressed as Based on Eq. (S6b), at T M we can define a fitting parameter ν as ν = T M D cf /D pt = T M C cf /C pt (S18) By using the four Eqs. (S16)-(S18) at T M , we can determine D pt , D cf , κ ρ , κ η and thus eight model parameters according to the dependency. Then according to Eqs. (S9) and (S15), the temperature-dependent γ sf (T ) and γ gb (T ) can be obtained. Normally we can fit γ sf (T ) and γ gb (T ) to make them close to the experimental data by giving λ TM gb and tuning ν. However, it is too thin for the real grain boundary width (usually several nanometers) when compare to the size of the particles/grains in the scale of micrometer, requiring much finer mesh in finite element simulations [14]. To deal with such a problem, a common solution is to manually set the diffusion width λ TM gb close to the scale of the particles/grains, e.g. 1 ∼ 2 µm for the particles/grains with an average size of several tens micrometer, while keep the γ gb identical. According to Eqs. (S16) and (S17), the barrier height of the potential part of the free energy density f pt = 12D (1 − η) 2 η 2 , which is proportional to the D, will decrease when wider γ gb is utilized to keep γ gb identical. In this case, the heat part of the free energy density f ht in Eqs. (S9) and (S15) should be modified as ξf ht by introducing a coefficient ξ to have the same scale as the potential part f pt . Fig. S5a shows the fitted γ sf and γ gb by setting λ TM gb = 2 µm (T M = 1700 K), ξ = 0.3 × 10 −3 , and ν = 1.01, as well as the experimental γ sf and γ gb .