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

TL;DR

This paper introduces a comprehensive 3D non-isothermal phase-field model for simulating microstructure evolution during selective laser sintering, capturing complex phenomena like melting, pore formation, and grain migration, validated through stainless steel 316L experiments.

Contribution

It presents a novel thermodynamically consistent non-isothermal phase-field model with efficient grain tracking for SLS, enabling detailed microstructure prediction and analysis.

Findings

Identified effects of laser power and scanning speed on microstructure.

Validated porosity evolution follows first-order kinetics.

Linked densification to specific energy input during SLS.

Abstract

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.

Equations

$$\begin{array}{rl}f(T,\rho ,\{{\eta }_{\alpha }\})=& \xi f_{\text{ht}}(T)\left(\underline{A}\rho +\underline{B}\sum _{\alpha }{\eta }_{\alpha }\right)+\underline{C}\left[{\rho }^2(1-\rho {)}^2\right]+\\ & \underline{D}\left[{\rho }^2+6(1-\rho )\sum _{\alpha }{\eta }_{\alpha }^2-4(2-\rho )\sum _{\alpha }{\eta }_{\alpha }^3+3{\left(\sum _{\alpha }{\eta }_{\alpha }^2\right)}^2\right]+\\ & \frac{1}{2}T{\kappa }_{\rho }{\left|\nabla \rho \right|}^2+\frac{1}{2}T{\kappa }_{\eta }\sum _{\alpha }{\left|\nabla {\eta }_{\alpha }\right|}^2,\end{array}$$

$$\begin{array}{rl}& f_{\text{ht}}(T)=T\int e_{\text{ht}}(T)\text{d}\left(\frac{1}{T}\right),\\ & \underline{C}={\underline{C}}_{\text{pt}}-{\underline{C}}_{\text{cf}}(T-T_{\text{M}}),\\ & \underline{D}={\underline{D}}_{\text{pt}}-{\underline{D}}_{\text{cf}}(T-T_{\text{M}}).\end{array}$$

$$e_{\text{ht}}(T)=(c_{\text{bk}}-c_{\text{at}})(T-T_{\text{M}})+\mathcal{L}\Phi (\tau ).$$

$$\dot{\rho }+\nabla \cdot {\mathbf{j}}_{\text{diff}}+\nabla \cdot {\mathbf{j}}_{\text{melt}}=0,$$

$$\dot{\rho }(\mathbf{r},t)=\nabla \cdot \left[\mathbf{M}\cdot \nabla \left(\frac{\partial f}{\partial \rho }-T{\kappa }_{\rho }{\nabla }^2\rho \right)\right].$$

$$\mathbf{M}={\left[{\mathbf{M}}_{\text{bk}}+{\mathbf{M}}_{\text{at}}+{\mathbf{M}}_{\text{sf}}+{\mathbf{M}}_{\text{gb}}\right]}_{\text{diff}}+\Phi (\tau ){\mathbf{M}}_{\text{melt}}.$$

$${\dot{\eta }}_{\alpha }(\mathbf{r},t)=-\mathbf{L}\left(\frac{\partial f}{\partial {\eta }_{\alpha }}-T{\kappa }_{\eta }{\nabla }^2{\eta }_{\alpha }\right).$$

$$\frac{\partial e}{\partial T}\left[\dot{T}(\mathbf{r},t)-\mathbf{v}\cdot \nabla T\right]+\frac{\partial e}{\partial \rho }\dot{\rho }(\mathbf{r},t)+\sum _{\alpha }\frac{\partial e}{\partial {\eta }_{\alpha }}{\dot{\eta }}_{\alpha }(\mathbf{r},t)=\nabla \cdot \left(\mathbf{k}\cdot \nabla T\right)+q(\mathbf{r}),$$

$$\mathbf{k}={\mathbf{k}}_{\text{bk}}+{\mathbf{k}}_{\text{at}}+{\mathbf{k}}_{\text{sf}}+{\mathbf{k}}_{\text{gb}}.$$

$$q(\mathbf{r})=-\beta P_0{p}_{xy}(x,y)p_z(\alpha ,\lambda ,z),$$

$$⟨\beta ⟩=\frac{3}{2}\frac{1-\epsilon }{\epsilon }⟨\frac{1}{d}⟩,⟨\lambda ⟩=\frac{3}{2}\frac{1-\epsilon }{\epsilon }⟨\frac{H_{\text{pb}}}{d}⟩.$$

$$-(\mathbf{k}\cdot \nabla T)\cdot \mathbf{n}=-\frac{k_{\text{sub}}}{H_{\text{sb}}}({T|}_{{\Gamma }_{\text{sub}}}-T_{\text{P}}),$$

$$M_{\text{path}}^{\text{eff}}=\frac{D_{\text{path}}^{\text{eff}}}{2(\underline{C}+\underline{D})},L^{\text{eff}}=\frac{G_{\text{gb}}^{\text{eff}}{\gamma }_{\text{gb}}}{T{\kappa }_{\eta }}.$$

$$V=\int \rho \delta ({\Omega }_i)\text{d}\Omega ,T_{\text{avg}}=\frac{\int \rho T\delta ({\Omega }_i)\text{d}\Omega }{\int \rho \delta ({\Omega }_i)\text{d}\Omega },$$

$$\dot{\epsilon }=-k\epsilon ,$$

$$\Psi =\frac{P}{hwv},\Omega =\frac{{\epsilon }_0-\epsilon }{{\epsilon }_0-{\epsilon }_{\text{min}}}.$$

Full text

1 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 [8, 5, 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 [10, 1, 11, 5]. 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 even sort of empirical knowledge [5, 19, 20]. There are already many computational works performed regarding the underlying physics during SLS, where powder-laser interaction [7, 21, 22, 23] and heat transfer [24, 25, 26, 27] have been massively investigated. Apart from those, the microstructure evolution has also gained great interests since it is directly related to mechanical properties of the final products like tensile strength, ductility, and fracture toughness. Up to now, however, it still remains a great challenge to simulate the microstructure evolution if the influences from all aspects are considered. During scanning, there is drastic difference in the thermal conditions among particles due to different exposure to the laser beam. Some of particles may even partially melted during the process [5, 28]. Therefore, binding mechanism for certain particle/grain may vary from the solid-state sintering to the liquid-state sintering, and even melting-solidification according to the intensity of its partial melting [8, 28, 29, 30]. For the same reason, very high temperature gradients and cooling rates also exist [31, 7], which make the mechanism of grain coalescence and coarsening during SLS deviate from that in conventional isothermal sintering.

In recent decades, the phase-field method has been utilized to simulate the microstructure evolution during sintering-related techniques because of its ability to model the complex spatial geometries without explicitly tracking the position of the surface and interface. One of the earliest attempts to combine such efforts into a phase-field sintering model is carried out by Wang et al., who utilized a conserved density field and a set of non-conserved orientation fields [32, 33]. Based on this model, some featured details, such as the dihedral angle, neck growth, shrinkage and the grain growth kinetics with porosity, have been captured and validated by experiments [34, 35, 36]. In recent research, Zhang et al. adopted this isothermal model to investigate the necking among several powders during the SLS [37], whereas the information of non-isothermal profile, i.e. inhomogeneous temperature field with high temperature gradient, is missing. Lu et al. developed a two-dimensional (2D) phase-field model to study the molten pool geometry, porosity, and grain structure during powder bed-based additive manufacturing [38]. This model presented a good coherency of porosity and grain structure with the experimental data. However, it is not accessible to the features like the surface morphology and grain geometry due to the limitation of 2D presentation. In brief, a three-dimensional (3D) non-isothermal phase-field model is still in demand to improve the understanding of the microstructure evolution during the SLS process.

In this work, we present a 3D non-isothermal phase-field model to investigate the microstructure evolution during the SLS process. The model is derived in a thermodynamically consistent way according to our latest work in [39] and numerically implemented with the finite element method (FEM). Besides the diffusion and grain boundary migration mechanisms, the role of the partial melting of the powders is also considered. A scenario with a conserved order parameter (OP) as well as a series of non-conserved ones is utilized to represent the microstructure evolution during SLS. The discrete element method (DEM) is used to generate the initial powder bed. The Welsh-Powell algorithm for the minimum coloring problem (MCP) is also applied to optimize the profile of OPs to reduce the computation cost due to large amount of variables. This model is demonstrated to be capable of capturing interesting phenomena which are not accessible to the conventional isothermal model. We further perform phase-field simulations on the SLS processing of stainless steel 316L powder, in which the temperature-dependent model parameters are readily extracted from the experimental measurement of surface and interface energies. The influences from the laser power and scanning speed on key features, such as the porosity, surface morphology, temperature profile, grain geometry and densification, are also discussed. It is hoped that the present work could enrich the modeling and computational toolkit which is practicable for the simulation of SLS-based additive manufacturing.

2 Modeling

2.1 Non-isothermal Phase-field Scenario

A segment of a powder bed is considered, as shown in Fig. 1a. A domain with the dimension of $250\times 500\times 100~{}$\mathrm{\SIUnitSymbolMicro m}$^{3}$ is utilized to simulate a powder bed. Both conserved and non-conserved order parameters (OPs) $\rho$ and $\{\eta_{\alpha}\}$ are employed to represent the powder bed with multiple particles. The conserved OP $\rho$ indicates the substance while the non-conserved OPs $\eta_{\alpha}$ distinguish particles with different crystallographic orientations [40, 34, 39]. $\rho=1$ and $\rho=0$ represent the substance and atmosphere/pore region, respectively. In each grain within the substance, only one of $\eta_{\alpha}$ takes unity and others are zero (Fig. 1b). These grains and grain boundaries have $\rho=1$ (we assume the density variation across the grain boundary is negligible). When $\rho=0$, no grain is present. This profile of OPs leads to the constraint $(1-\rho)+\sum_{\alpha}\eta_{\alpha}=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,\rho,\{\eta_{\alpha}\})$, derived through the Legendre transform of the internal energy and entropy [41], is formulated as

[TABLE]

where

[TABLE]

The eight model parameters ($\underline{A},\underline{B},\underline{C}_{\text{pt}},\underline{C}_{\text{cf}},\underline{D}_{\text{pt}},\underline{D}_{\text{cf}},~{}\kappa_{\eta}$ and $\kappa_{\rho}$) in Eq. (1) are obtained from the experimentally measured temperature-dependent surface and grain-boundary energies ($\gamma_{\text{sf}}$ and $\gamma_{\text{gb}}$) and the width of the grain boundary. Fitting coefficient $\xi$ is utilized to help fitting with the experimental result (see Supplementary Information). $e_{\text{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_{\text{M}}$ of the system as zero)

[TABLE]

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 $\mathscr{L}\Phi(\tau)$, where the sufficiently smooth interpolation function $\Phi(\tau)$ takes one when $\tau=T/T_{\text{M}}\rightarrow 1$ and to zero when $\tau\rightarrow 0$.

The temporal evolution of $\rho$ is derived from the mass conservation, i.e.

[TABLE]

where ${{\mathbf{j}}_{\text{diff}}}$ and ${{\mathbf{j}}_{\text{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_{\text{sf}}$ and energy $\gamma_{\text{sf}}$, i.e. ${{\mathbf{j}}_{\text{melt}}}\sim-\Phi(\tau)\mathbf{M}_{\text{melt}}\cdot\nabla(A_{\text{sf}}\gamma_{\text{sf}})$. According to Young-Laplace equation in which the on-site chemical potential is proportional to $A_{\text{sf}}\gamma_{\text{sf}}$, we can thereby formulate ${{\mathbf{j}}_{\text{melt}}}$ into a diffusion-like form. ${{\mathbf{j}}_{\text{diff}}}$ can be obtained directly from Fick’s law. In this way, Eq. (3) can be rearranged into the form of the Cahn-Hilliard equation

[TABLE]

Here the mobility tensor $\mathbf{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) [40, 42, 43, 34], but also from the gain due to the partial melting, i.e.

[TABLE]

In this regard the contribution of the partial melting is treated more like an enhanced surface diffusion when $\tau\rightarrow 1$. To represent such contribution more properly, a ${{\mathbf{j}}_{\text{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 $\{\eta_{\alpha}\}$ is governed by Allen-Cahn equation with the corresponding mobility tensor $\mathbf{L}$, i.e.

[TABLE]

The laser scanning process is modeled by moving the Gaussian distributed heat source on the powder-bed surface at a speed of $\mathbf{v}$. The kinematic equation of temperature is formulated as

[TABLE]

where the thermal conductivity $\mathbf{k}$ adopts spatial distribution as $\mathbf{M}$ in Eq. (5) i.e.

[TABLE]

$q(\mathbf{r})$ represents the volumetric energy deposition due to the radiative energy flux of the laser and can be formulated as [7, 6]

[TABLE]

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}(\alpha,\lambda,z)$ is a penetration function, which takes the form proposed by Gusarov et al. [6]. The hemispherical reflectivity $\alpha$ indicates the influence of the powder material. The attenuation coefficient $\beta$ and the optical thickness $\lambda$ reflect the influence from the powder bed structure. For a loosely packed powder bed with particles of distributed diameters $d$, a thickness of $H_{\text{pb}}$, and a porosity of $\varepsilon$, the effective value of these two parameters are calculated as

[TABLE]

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.

[TABLE]

where $\mathbf{n}$ is the normal vector of the boundary $\Gamma_{\text{sub}}$. $H_{\text{sb}}$ and $H_{\text{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_{\text{P}}$ as the pre-heating temperature.

2.2 Optimized profile of order parameters

It is anticipated that a general powder bed would have a spatial and temporal conserved OP $\rho(\mathbf{r},t)$ in combination with the multiple non-conserved OPs $\{\eta_{\alpha}\}=\{\eta_{1}(\mathbf{r},t),\eta_{2}(\mathbf{r},t),\dots\eta_{N}(\mathbf{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 next-nearest neighbor in hexagonal close packing (inset of in Fig. 1c) is considered as the criterion of adjacency. Then, the map $G(\mathbb{V},\mathbb{E})$ is established based on the adjacency matrix, which consists the unique particle $\nu_{i}(\mathbf{r},R)$ as the vertices $\mathbb{V}\{\nu_{i}\}$ and corresponding degree of adjacency $\mathbb{E}\{\epsilon_{i}\}$, i.e. there are $\epsilon_{i}$ particles are adjacent to the one indexed as $\nu_{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 $\mathrm{\SIUnitSymbolMicro 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.

3 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 $\mathbf{L}$ and $\mathbf{M}$ are formulated as scalars and adopted from the self-diffusivity $D_{\text{path}}^{\text{eff}}$ through possible paths (path $=$ bk, vp, sf and gb) and grain boundary mobility $G_{\text{gb}}^{\text{eff}}$ [50, 39, 36], i.e.

[TABLE]

Details of the material-related quantities are elaborated in Note 1 (Supplementary Information), and the procedure to obtain model parameters ($\underline{A},\underline{B},\underline{C}_{\text{pt}},\underline{C}_{\text{cf}},\underline{D}_{\text{pt}},\underline{D}_{\text{cf}},~{}\kappa_{\eta}$ and $\kappa_{\rho}$) are elaborated in Note 2 (Supplementary Information). The thickness of the powder bed with SS316L particles and the SS316L substrate is $H_{\text{pb}}=65~{}$\mathrm{\SIUnitSymbolMicro m} and $H_{\text{sb}}=1000~{}$\mathrm{\SIUnitSymbolMicro m}, respectively. The SLS atmosphere is argon. The powder bed consists of particles with distributed diameters $d$ ranging from $20$ to $50~{}$\mathrm{\SIUnitSymbolMicro m}. The distribution of $d$ is shown in Fig. S1b. (Supplementary Information). According to Eq. ([10](#S2.E10)), the effective attenuation coefficient and optical thickness of the powder bed are calculated as $\langle\beta\rangle=~{}0.089$ and $\langle\lambda\rangle=5.8$. The Gaussian distribution of the laser beam $p_{xy}(x,y)$ has the full width at half maximum $\text{FWHM}=100~{}$\mathrm{\SIUnitSymbolMicro m}. Note that the surface integral of $p_{xy}(x,y)$ reaches $75.8\%$ of $P_{0}$ in the region with a diameter of FWHM, and reaches $99.5\%$ in the region with a diameter of $2\times\text{FWHM}$. So we regard the region $2\times\text{FWHM}$ as the nominal beam spot. Temperature of the powder bed is initialized as $T/T_{\text{M}}=0.4$ with the melting point $T_{\text{M}}=1700~{}$\mathrm{K}$$.

Fig. 2a presents the microstructure evolution obtained at a laser power of 20 W and a scanning speed of 100 $\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{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_{\text{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_{\text{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_{\text{M}}$. In the region with a maximum temperature lower than $T_{\text{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~{}$\mathrm{K}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{m} and $100~{}$\mathrm{K}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{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.3b4). 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. 3b4-3b7; 3b8-3b10; and 3b3, 3b12-3b11). It is similar for the case when scanning speed is fixed and the laser power is increased (comparing Fig. 3b1-3b4; 3b5, 3b8, 3b12; and 3b6, 3b10-3b11).

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.

4 Discussion

4.1 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 [24, 6, 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 $\rho$ is introduced to assign different thermal properties to the substance and pores separately. The microstructure descriptor $\rho$ and $\{\eta_{\alpha}\}$ 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\times\text{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 $\mathrm{\SIUnitSymbolMicro 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~{}$\mathrm{W} and $v=100~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}, we can only see some particles are partially melted. When increasing the laser power (e.g. $P_{0}=30~{}$\mathrm{W} and $v=100~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s} in Fig. 4b) or decreasing the scanning speed (e.g. $P_{0}=20~{}$\mathrm{W} and $v=80~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{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~{}$\mathrm{W} and $v=100~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s} in Fig. 4b) or increasing the scanning speed (e.g. $P_{0}=20~{}$\mathrm{W} and $v=150~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{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_{\text{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\times\text{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.

4.2 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~{}$\mathrm{K}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{m} within the particles and $100~{}$\mathrm{K}\mathrm{/}\mathrm{\SIUnitSymbolMicro}\mathrm{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~{}$\mathrm{W} and $v=100~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}. Here the volume and the average temperature of the $i\text{-th}$ particle/grain are calculated as

[TABLE]

where $\delta(\Omega_{i})$ is the delta function which takes unity inside the domain $\Omega_{i}$ of the $i\text{-th}$ particle/grain. $\Omega$ represents the simulation domain. The volume variation $\Delta 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_{\text{M}}$. After reaching the maximum $T_{\text{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 $\mathrm{\SIUnitSymbolMicro s}$), they are in a polyhedral shape with arched surfaces (g${}_{\text{A}}$ ,g${}_{\text{D}}$ and g${}_{\text{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${}_{\text{B}}$). Grain C is located closest to the center of the partial melting region, whose $T_{\text{avg}}$ firstly rises above $T_{\text{M}}$ then quickly drops. It eventually turns into a flat grain (g${}_{\text{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_{\text{avg}}$. In the end, it also turns into a polyhedral shape (g${}_{\text{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~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}, the grain is polyhedral at $P_{0}=20~{}$\mathrm{W} (g3), but is round with a top-surface necking (g5) at $P_{0}=10~{}$\mathrm{W}. For the cases of $P_{0}=30~{}$\mathrm{W}, $v=100~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s} and $P_{0}=20~{}$\mathrm{W}, $v=80~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}, flat grains with a smooth top and a polyhedral bottom are obtained (g1 and g2). We can also observe the fluctuation of the volume variation in the case of $P_{0}=30~{}$\mathrm{W}, $v=100~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}$$ when $T_{\text{avg}}$ is approaching its maximum, which is close to the melting point. This may be due to the contribution from the partial melting term $\Phi(\tau)\mathbf{M}_{\text{melt}}$.

4.3 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

[TABLE]

where $k$ is defined as the sintering rate constant. In Fig. 6a we present the temporal evolution of $\ln({\varepsilon/\varepsilon_{0}})$ of the three selected segments and the whole powder bed with an initial porosity $\varepsilon_{0}$ under the processing parameters of $P_{0}=20~{}$\mathrm{W}, $v=100~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}. The average temperature $T_{\text{avg}}$ of each segment is also calculated according to Eq. (13). These three segments are with the same geometry (length 100 $\mathrm{\SIUnitSymbolMicro m}$, width 250 $\mathrm{\SIUnitSymbolMicro 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_{\text{avg}}$ rises and $\ln({\varepsilon/\varepsilon_{0}})$ is almost linearly decreased there. Once the laser beam is about to leave the certain segment, the $T_{\text{avg}}$ reaches its maximum then start to drop, while $\ln({\varepsilon/\varepsilon_{0}})$ is slowly decreased. In contrast, $\ln({\varepsilon/\varepsilon_{0}})$ of the whole powder bed presents an approximately linear trend, with a fitted slope as $k=7\times 10^{-5}~{}$\mathrm{s}^{-1}, as shown in Fig. [6](#S4.F6)a. This demonstrates that Eq. ([13](#S4.E13)) 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](#bib.bib19)]. According to Fig. [6](#S4.F6)b, we can see the final microstructure ($t=5000~{}$\mathrm{\SIUnitSymbolMicro 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({\varepsilon/\varepsilon_{0}})$ of the whole powder bed under different processing parameters. Segment II is selected to present the influence of processing parameters on temporal $T_{\text{avg}}$, which shows a similar trend as in Fig. 6a. The approximately linear trend of $\ln({\varepsilon/\varepsilon_{0}})$ is still obvious. Both increasing the laser power and decreasing the scanning speed result in the increase of maximum $T_{\text{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~{}$\mathrm{W}, $v=100~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s} as the reference and compare it with cases with other processing parameters. For the cases with fixed $v=100~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s} and $P_{0}=10$, $15$ and $20~{}$\mathrm{W}, densification rarely occurs and only necking with tunnel-like open pores form, also manifesting the characteristics of the initial stage of sintering. For the case of $P_{0}=30~{}$\mathrm{W} and $v=100~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}, however, the microstructure is significantly densified with packed grains and almost no pores, showing the characteristics of the final stage of the sintering where the grain coarsening dominates. For the case of $P_{0}=20~{}$\mathrm{W} and $v=80~{}$\mathrm{m}\mathrm{m}\mathrm{/}\mathrm{s}, the microstructure densification occurs, but there are still some round pores near the segment margin, indicating the transition from the intermediate stage to the final stage of sintering.

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 $\Psi$ and densification factor $\Omega$ as

[TABLE]

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. $\varepsilon_{\text{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 $\text{FWHM}=100~{}$\mathrm{\SIUnitSymbolMicro 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_{\text{pb}}=65~{}$\mathrm{\SIUnitSymbolMicro m} and $\varepsilon_{\text{min}}=3\%$. As shown in Fig. 7a, the relation $k$ vs $\Psi$ can be linearly fitted with a slope of $9.51\times 10^{-3}$. Similarly, linear fitting of $-\ln(1-\Omega)$ vs $\Psi$ 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 $\mathrm{\SIUnitSymbolMicro 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 $\Omega$ in Fig. 7c. Following the isoline of the specific energy input $\Psi$ (dot-dash lines in Fig. 7c), increasing laser power/scan speed actually leads to the increase of the densification $\Omega$. And when higher $\Psi$ are kept, the faster increase of $\Omega$ 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.

5 Method

5.1 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 $\mathbf{r}$ and radius $R$ of each unique particle is recorded as a vertex $\nu_{i}(\mathbf{r},R)$, and added into the vertices list $\mathbb{V}\{\nu_{i}\}$ for further optimization.

5.2 Implementation of finite element method

The model is numerically implemented within the MOOSE framework by finite element method (FEM) [56]. 8-node 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 $\mu$ [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.

5.3 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.

5.4 Data availability

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

6 Acknowledgement

The support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 743116) and the Profile Area From Material to Product Innovation – PMP (TU Darmstadt) is acknowledged. The authors also greatly appreciate their access to the Lichtenberg High Performance Computer and the technique supports from the HHLR, Technische Universität Darmstadt.

7 Author Contributions

Y.Y. performed the phase-field simulations, data processing, and manuscript writing. O.R. performed the powder bed generation and contributed to the manuscript writing. Y.B. supported numerical calculations in HPC. B.-X.X. supervised the project, analyzed the results, and revised the manuscript. M.Y. contributed by intensive consultation and discussion. All authors reviewed and approved the manuscript.