Microstructure and Elements Concentration of Inconel 713LC during Laser Powder Bed Fusion through a Modiﬁed Cellular Automaton Model

: In this study, a hybrid ﬁnite element (FE) and cellular automaton (CA) model is developed to explore crystallization behavior and alloying of Inconel713LC during Laser powder bed fusion. A cellular automaton model is considering the surface nucleation, equiaxed bulk nucleation, and grain growth kinetics. In addition, the equation for solute diffusion is coupled with a cellular automaton model to simulate the IN713LC elements segregation. During the phase change, the non-equilibrium segregation model is applied to insert the effect of ultra-fast solidiﬁcation happening during LPBF. It is found that, during LPBF processing of IN713LC, the micro segregation of Nb, Ti, and C is accrued at the grain boundaries. It is further shown that the micro segregation intensity depends on the solidiﬁcation speed, which is determined in turn by the laser heat input. In particular, a lower laser heat input increases the solidiﬁcation speed and results in a more uniform solid phase, thereby reducing the risk of crack formation. Finally, using a comparison between simulation results and experimental observation, it was shown that the proposed model successfully predicts the bulk element concentration of IN713LC after laser melting.


Introduction
IN713LC, a nickel alloy with a composition of Cr-Al-Mo-Ti-Nb-Zr, is known for its good fatigue resistance, excellent mechanical properties, superior oxidation resistance, and enhanced resistance to degradation under harsh operating conditions [1][2][3][4]. Although the Ni-based superalloys faced many challenges with the LPBF process, still many studies have shown that most of them such as IN718 [5], IN625 [6], and Hastelloy X [7] can be successfully fabricated by laser powder bed fusion (LPBF). However, due to its crackprone nature, the LPBF processing of IN713LC is extremely challenging [8,9]. Due to the rapid rate of solidification during LPBF, experimental investigations provide only limited insights about LPBF of IN713LC. Thus, in exploring the mechanisms which render the LPBF processing of IN713LC so challenging, simulation models are widely preferred. Cellular Automaton (CA) is recognized as one of the best simulation methods for examining ultrafast solidification processes such as LPBF [10][11][12][13]. Ao et al. [14] used CA to simulate the microstructure of AlSi10Mg with a moving molten pool by selective laser melting. They have used the thermodynamic driven growth kinetics of the solid/liquid and composition. In addition, the growth rate was computed by the Kurz-Giovanola-Trivedi model [15] and they showed that the cooling rate is around 105-106 K/s. In addition, they proved that, during the solidification, the equiaxed grain formation increases with increasing the pre-heating temperature or reducing the scanning speed. A two-dimensional CA model is developed by Zinoviev et al. [16] to evaluate grain structure of alloys during SLM. They proved that the CA modelling can provide an accurate prediction of the alloy microstructure made by SLM. In addition, the microstructure of IN718 is predicted by a

Finite Element Model
The transition governing equations used in this study are expressed as follows [21]: where ρ is density, h is the enthalpy, → V is the velocity field, T is temperature, . Q L is the laser heat source. In the momentum equation, p denotes pressure, µ is viscosity, and → V ALE is the mesh velocity field from the Arbitrary Lagrangian-Eulerian (ALE) model [22]. In addition, S b represents the buoyancy forces and S v is the body force that is expressed by Carman-Kozeny approximation obtained from Darcy's law as [23][24][25]: where β is the thermal expansion coefficient, → g is the gravity vector, Tm is the IN713LC melting temperature, f l is the liquid fraction, ε is a small number to prevent division by zero during the solidification, and C is a constant number which is set here to 106. The laser heat source is based on volumetric laser-powder and laser-bulk matter (melted or solidified IN713LC) interactions and defined as [24]: Here, the liquid fraction, f l , and the laser-melt interaction, . q m (r) are expressed by [24]: In addition, . q p (r) is the laser-powder interaction which is modified by Beer-Lambert law [26] and has the following form: where γ is the angle between laser direction and normal melt surface, r 0 is the laser radius, P l is the laser power, R is the IN713LC reflection coefficient, r is the distance between any point and laser center incident point, α is in-depth absorption coefficient, and r * 0 is the effective laser interaction radius where, due to radiation between powders, r * 0 > r 0 . The effect of latent heat during phase change is included in specific heat expression as [24,27,28]: where ∆H m is the latent heat of fusion and ω is the smoothed function to express the fraction of phase during the phase change and can be expressed as [28]: Thermal conductivity, k, and density, ρ, of bulk IN713LC are respectively defined as [28]: where k s and k l are the thermal conductivity of IN713LC at solid and liquid phases. In order to model the powder physical properties, Sih [29] effective conductivity is used as: where ϕ is the porosity of the powder layer, k g is the thermal conductivity of atmosphere gas, and k r is the radiation thermal conductivity between individual powder particles, which is defined as [29]: where σ is the Stefan-Boltzmann constant, and D p is the average powder size. To find powder surface emissivity, powder is considered as a porous media including particles and cavities. In this condition, the powder media emissivity is obtained from summation of emissivity of particles and cavities located at the powder surface as [29]: where ε H is the emissivity of surface cavities and A H is the fraction of cavities on the powder surfaces. These two parameters can be defined as [29]: All other powder's physical properties, ∅ p , can be expressed as a function of corresponding IN713LC property, ∅, and the porosity of the powder layer, ϕ, as [30]: Both Marangoni and surface tension forced are imposed as boundary conditions on melt surface as: where n is the unit vector normal to the melt surface and γ is the surface tension coefficient. An arbitrary Lagrangian Eulerian (ALE) method [22] is implemented in order to simulate the free melt surface evaluation. The moving mesh velocity can be expressed as: where v sh is the powder volume shrinkage speed during the melting and can be expressed as:

Microstructure Model
In the CA model, the number of nucleation sites formed during the solidification is conventionally assumed to have a Gaussian distribution as a function of the undercooling temperature and can be modeled using the following equation [10]: where ∆T is the undercooling temperature, ∆T 0 max is the maximum nucleation undercooling temperature, ∆T σ is the standard deviation of the undercooling temperature, and n max is the maximum nucleation density. The total density of the nucleation sites can thus be modeled as Since the IN713LC alloy system is dominated by Ni (~75 wt.%), the undercooling temperature, ∆T, can be assumed to consist of only thermodynamic undercooling ∆T T . That is, The thermodynamic undercooling, ∆T T , can be expressed as [10]: where c p , ∆H, P t , and Iv are the specific heat, latent heat, Peclet number, and Ivantsov function, respectively. The Ivantsov function has for [10]: In addition, the Peclet number is defined as [10]: where v growth is the liquid/solid growth speed, α is the thermal diffusivity, and R is the dendrite tip radius. According to Trivedi [31] and Kurz and Fisher [32], the growth speed and dendrite tip radius are related via the following stability criterion: where Γ is the Gibbs-Thomson coefficient. From Equations (28) and (29), the growth rate can be obtained as v growth In modeling the element segregation behavior during solidification, the governing equation for solute diffusion is given as [33]: where D i is the solute diffusion coefficient of component i (Cr, Mo, Nb, Al, Ti, and C in the present case). In addition, the third term in Equation (31) is the solute source term which is segregated at the solid-liquid interface and C * L,i and C * S,i are the solute concentrations of component i at the solid-liquid interface, respectively. Under equilibrium conditions, the two solute concentration terms at the solid-liquid interface are related as follows: where k i is the equilibrium partition coefficient of element i. However, for the rapid solidification rate considered in the present study, the non-equilibrium relation between the element concentrations in the liquid and solid phases is given as where a 0 is the length of atomic dimensions of the components and R s is the solidification rate, as computed by the CA model. From Equation (33), when the solidification speed is slow, R s → 0 and C * s,i = k i C * L,i . However, when the solidification speed is very fast, R s → ∞ and C * s,i = C * L,i . In other words, no element segregation occurs since the elements are captured by the solid phase before they can diffuse to the liquid phase.
In the CA model, the simulation domain is divided into cells with a uniform size of 0.2 µm × 0.2 µm. Each cell is characterized by two variables, namely the state (i.e., solid or liquid) and the crystallography orientation. In initializing the model, the system is considered to be at the solidus temperature, and hence all of the cell states are set as 1. During the simulation, if temperature of the CA cell passes from liquid temperature, its state and crystallography are changed to 0. Notably, the cell state may change from 0 or 1 under the following two conditions:

•
Nucleation occurs in the cell; • The cell is captured by a solid cell.
The density of new nuclei formed during the solidification process is calculated by Equation (23) for simulation domain at boundaries or inside the melt. The nucleation densities are then multiplied by the total number of cells with a state of 1 in order to determine the number of new nucleation sites. That is, where N as is the total number of cells with state 1 at the system boundary, and N av is the total number of cells with state 1 inside the bulk liquid. Furthermore, ∆T 0 vmax is considered to be 5 • C lower than ∆T 0 smax . In performing the simulations, a random number, ri, with a value between 0 and 1 is assigned to each CA cell and nucleation is assumed to occur if the following condition is then satisfied: Following nucleation, a second random process is used to generate and assign a crystallography orientation number, q = [0−64], to the corresponding cell. Note that a large number of possible crystallography orientations is deliberately considered here in order to avoid impinging grains with the same orientation during the solidification process. In addition, the growth length of cell i with regard to its neighbor j at time t is calculated as Liquid neighbor capturing is also regarded as a random process. In particular, neighbor j is captured by cell i if the following condition is satisfied: where ∆x is the cell length and θ is the crystallographic orientation, which is calculated from the state of cell i as follows: Equations (37) and (38) show that the grain growth speed is a function of the grain orientation. In particular, the grain orientation affects the capturing rate and then drives the grain growth phenomena. After cell i captures cell j, the orientation of cell j is assigned to that of cell i. The solid fraction of cells, f S , is then calculated at each time step using the following equation: Considering the need for convergence and numerical stability, the time step in the simulation process is limited as follows: The 3D-FE model including all thermal equations is explicitly implemented in COM-SOL Multiphysics. A microstructure model is also implemented in the in-house written code using FORTRAN. Heat loss radiation and convection boundary conditions are applied for the surface of the system. The Laplace smoothing algorithm is used to control the displacements of nodes. During the ALE model, if meshes deform too much, the map from mesh coordinates to spatial coordinates may enhance the ill-condition. In this condition, the auto re-mesh is applied. The powder layer thickness is considered 40 µm where the powder porosity, ϕ, is assumed to be 50%. For all the simulations, after laser track end up, simulation continues until the system became solid. In implementing the CA model, the IN713LC alloy was considered to be a 7-element system consisting of 12.1 Cr, 4.10 Mo, 2.0 Nb, 6.2 Al, 0.77 Ti, 0.05 C, and balance Ni. Note that the C component of the alloy was retained in the model even though its concentration is extremely low due its importance in increasing the formation of hot cracking. Furthermore, the diffusion coefficients in the solid phase are considered five orders smaller than liquid phase and no-flux boundary conditions were applied at all the boundaries. Table 1 summarizes the parameters applied in the micro-segregation model for the different elements of the IN713LC alloy. The physical parameters employed in the CA model are listed in Table 2.

Laser Processing
The printing trials were performed on a Tongtai Taiwan, Kaohsiung city, AM-250 LPBF machine equipped with a Nd-YAG laser with a wavelength of 1064 nm, a focal beam diameter of 100 µm, and a Gaussian irradiation profile. The laser provided a maximum power of 500 W and a maximum scanning speed of 2000 mm/s. A single track experiment was conducted on a substrate made up of IN713LC. A total of nine single tracks were deposited on the substrate for laser energy densities of 360, 280, and 210 J/m, respectively. A layer thickness of 40 µm was maintained for single-track deposits with powder. To capture the statistical variation, three identical single tracks for each parameter is conducted. The single track length was kept constant at 10 mm for all single track experiments.
The trials used a zig-zag scanning model with a 67 • rotation of the scanning direction between consecutive layers. To avoid oxidation during the LPBF process, the experiments were performed in a nitrogen atmosphere with an oxygen content of less than 1000 ppm. Based on a series of preliminary experiments, the scanning process was performed using three different laser energy densities of 360, 280, and 210 J/m, respectively. For each energy density, three cubic specimens were printed with dimensions of 10 × 10 × 10 mm 3 . Following the LPBF process, the samples were cut from the base plate with a wire and mounted on epoxy resin. The mounted samples were ground progressively with SiC sandpaper from a grit size of P240 to a final grit size of P3000, and were then further polished with a diamond suspension to a final particle size of around 0.3 µm. The polished samples were cleaned in ethanol solution and then dried. The microstructures and element compositions of the samples were examined using a scanning electron microscope (SEM, ZEISS Supra 55) equipped with energy dispersive X-ray spectroscopy (EDS) and electron backscattered diffraction (EBSD). In addition, the crack density was calculated using ImageJ software.

Results and Discussion
The 3D temperature distribution and surface pattern under laser power 250 W and scanning speed 700 mm/s, obtained from simulation, is shown in Figure 1a. In addition, Figure 1b,c compares the height map, surface morphology, and grain structure between simulation and experimental observations. The track width, melt penetration depth, and track bead height are obtained by simulation equal to 135 µm, 93 µm, and 56 µm, respectively. In addition, using experimental observation, the average track width, melt penetration depth, and track bead height are determined to be 143 µm, 83 µm, and 51 µm, respectively. Furthermore, grain structure shows the same trend between simulation and experimental results when the grains elongated toward melt pool center. Both experimental and simulation models show formation of equiaxed grains at the center of the melt pool. The average grain size is calculated using ImageJ software for experimental and simulation results equal to 56 µm and 49 µm. The comparison proved that the single-track geometry matches well with the experimental observations. The Q/v, the ratio of laser power to laser travelling speed, is known as the heat input [39]. Here, the microstructures and element compositions of the IN713LC samples fabricated using heat inputs of 360, 280, and 210 J/m were simulated using the present model. These heat inputs are obtained at laser powers 250 W, 200 W, and 150 W, while the laser scanning speed remained at~700 mm/s. Table 3 lists the results obtained from Equation (33) of the CA model for the average non-equilibrium partition coefficients of Cr, Mo, Nb, Al, Ti, and C during solidification given different values of the laser heat input. It is seen that the different elements have different solubilities in the IN713LC matrix. In particular, the C, Nb, Ti, Al, and Mo elements have relatively lower solubilities in the solid matrix, and hence tend to segregate in liquid. By contrast, Cr has a higher solubility, and thus tends to be trapped in the solid. It is noted that C has the highest diffusion coefficient of all the components in IN713LC and therefore moves through the solid-liquid interface easier than the other elements. The rapid diffusion of C atoms can be attributed to their small size, which allows them to sit interstitially within the Ni-super alloy lattice and hence move faster by jumping from one interstice to another. Comparing the equilibrium partition coefficients in Table 1 with the non-equilibrium element partition coefficients in Table 3, it is also seen that their values tend toward unity as the heat input reduces, i.e., the cooling rate increases. The average cooling rate during the solidification can be expressed by: (41) where N is the number of cells which experience phase change during the LPBF process, T i is the temperature of cell i, andt i is the total time that cell i remains in the liquid phase during the cooling. The average cooling rates during the solidification were calculated equal to 4.3 × 10 4 , 5.5 × 10 4 , and 7.3 × 10 4 K/s, for heat inputs 360, 280, and 210 J/m, respectively. In other words, a more uniform element concentration is obtained under lower energy densities. This finding is reasonable since, as the energy density decreases, the cooling rate and solidification speed increase, and hence the elements have less time to diffuse from the solid phase to the liquid phase. Figures 2-4 show the simulation results for the grain structure and Nb, Ti, and C concentrations after laser melting for laser heat inputs of 360, 280, and 210 J/m, respectively. Note that Nb, Ti, and C are deliberately selected here since they have the lowest partition coefficients among all the elements in IN713LC, as shown in Table 3. In general, the results reveal that a lower heat input leads to a more dominant columnar grain growth structure along the solidification direction. In addition, a higher heat input results in a greater number of equiaxed grains above the columnar grains (see Figure 2, for example, for the maximum heat input of 360 J/m). The greater volume fraction of equiaxed grains can be attributed to the lower temperature gradient produced at higher heat inputs (i.e., 4.5 × 10 5 K/m at an energy density of 360 J/m, compared to 5.8 × 10 5 and 7.8 × 10 5 K/m at heat input of 280 J/m and 210 J/m, respectively). More specifically, a lower temperature gradient results in the formation of equiaxed grains in front of the columnar grains, whereas a higher temperature gradient promotes directional growth.   The element distribution results in Figure 2 also show an obvious segregation of the Nb, Ti, and C elements at the grain boundaries. However, as shown in Figures 3 and 4, such micro-segregation effect is reducing under lower heat inputs of 280 J/m and 210 J/m, respectively. Results illustrated in Figure 4 show that, for low heat input 210 J/m that leads to a high cooling rate, there is almost no segregation during solidification. Hence, it can be inferred that LPBF processes performed under higher laser heat input (i.e., slower solidification speeds) will produce more intensive micro segregation of the elements. This intensive micro segregation may result in carbide or other secondary phase formation at the grain boundaries especially with the presence of Nb and Ti. Nb is founded as a strong carbide former [19]. In addition, based on the previous research, Ti is also known to control the formation of carbide during the solidification [20]. Then, segregation of both Nb and Ti during the solidification is important and even a small difference may result in different conditions.
A detailed inspection of the simulation results presented in Figures 2-4 shows that not all the grain boundaries exhibit the same amount of segregation. In general, grains with a larger size inevitably accrue a greater amount of segregation in front of their solid-liquid interface. When the segregated liquid is confined by other grains in the upstream region of the melt pool, the segregated elements are confined at the interface between them, and hence a region of high element concentration is formed. This is why most of the segregation happens at the boundary of large grains especially for heat inputs 360 J/m. Although the present model cannot simulate crack formation but is based on previous studies [18][19][20], there is a correlation between elements' segregations and crack formation, which shows that the present model can predict the crack formation qualitatively. According to the simulation results, the faster cooling rate suppresses element segregation at the grain boundaries and generates a more uniform element distribution. As a result, the formation of an undesired phase, which leads to crack initiation, is reduced. To verify the role of element segregation in prompting crack formation in the present samples, six samples in the form of single track or bulk cubic samples are made with different heat inputs which results in a different cooling rate, and then crack density is calculated for these three cases.  For the single track experiment, the crack density and average crack length reduce from 0.062% and 63 µm for the heat input 360 J/m to 0.006% and 15 µm for a heat input of 210 J/m. For bulk cubic samples, the crack density and average crack length reduce from 0.103% and 79 µm for the heat inputs of 360 J/m to 0.029% and 62 µm for an energy density of 210 J/m. In addition, a detailed EDS analysis was performed at selected points inside and outside of the cracks formed in the sample built with a high heat input of 360 J/m. The analysis points and corresponding EDS results are presented in Figure 5g,h. As predicted by the simulation results, both crack regions exhibit a high segregation rate of Nb, Ti, Al, and C. In particular, the Ti, Nb, Al, and C contents increase from 0.77, 2.0, 6.2, and 0.05 wt.% in the original IN713LC powder to 24.42, 11.54, 16.08, and 21.54 wt.%, respectively, within the crack. The high concentration of metallic elements implies the formation of MC carbides within the crack, and is hence consistent with the simulation results for the non-equilibrium element partition coefficients listed in Table 3.
In general, the simulation and experimental results presented above indicate that a lower heat input (i.e., a faster solidification rate) is beneficial in suppressing the segregation and then crack formation during the LPBF process of IN713LC components. Accordingly, a further sample was printed using a laser power of 120 W and a scanning speed of 700 mm/s, corresponding to a reduced energy density of around 170 J/m. Figure 6 presents an OM image of the built sample. It is seen that the surface contains neither cracks nor visible pores. In other words, the feasibility of LPBF for the processing of IN713LC is confirmed given a suitable low value of the laser energy density. It is noted, however, that reducing the energy density below a certain critical value may potentially induce new defects such as 'void formation' due to melt pool shrinkage under low energy conditions [40]. Here, to eliminate the sub-micrometric spherical pores at low energy density, the scan spacing was reduced to compensate for the shrinking of the melt pool. Accordingly, determining the optimal value of the energy density is an essential requirement for future studies. The same conclusion is obtained by Rashid et al. [41], where the heat input and energy per layer is found to be a critical parameter on the density and microstructure of final products during the LPBF process.

Conclusions
In conclusion, a modified CA model has been developed to predict the microstructure and micro segregation behavior of IN713LC nickel alloy during the solidification stage of the LPBF process. The below conclusions can be extracted from the results: (1) Solidification speed can be controlled by heat input. Faster cooling rate and solidification speed are obtained by decreasing the laser heat input. (2) The simulation results show that the micro segregation phenomenon is a diffusioncontrolled process, in which the elements with a lower partition coefficient and higher diffusivity experience a higher rate of segregation. (3) It is shown that element segregation is enhanced under a slower solidification rate since the elements spend a longer time in their respective precipitation windows and thus have sufficient time to move from the solid phase to the liquid phase. (4) As the laser heat input reduces, the solidification speed increases and the element partitioning coefficients approach unity. Consequently, a more uniform solid phase is formed with only minimal segregation at the grain boundaries. (5) It has been shown that the crack length and crack density decrease with a reducing laser heat input.
Notably, the experimental results carried out using the concepts of the present model have demonstrated the feasibility of fabricating crack-free IN713LC components using the LPBF process given a suitable low value of the energy density (e.g., 170 J/m in the present case). For future work, the present model can be used in machine learning and digital twining models and help to predict elements concentration and optimization of the LPBF process.