A multi-grid Cellular Automaton model for simulating dendrite growth and its application in additive manufacturing

Fig. 1a illustrates the multi-grid method in 2D for simplicity. First of all, similar to the conventional CA models, the calculation domain is discretized into uniform cubic cells, i.e. the coarse grid (hereafter called as parent grid), and the cell size is $d_{CELL}$. In the parent grid, each cell is associated with a set of variables ($F_{s}$, $C_{l}$, $g$, $S$): $F_{s}$ is the solid fraction of this cell, $C_{l}$ is the solute concentration in the liquid, $g$ is the identification number of the dendrite, to which the cell belongs, and $S$ is the status of this cell ($S=-1$: solidified cell; $S=0$: interfacial cell; $S=1$: liquid cell; $S=2$: potential interfacial cell). The interfacial cells (green cells in Fig. 1a) and potential interfacial cells (yellow cells in Fig. 1a) are tracked and further discretized into smaller uniform cubic cells (hereafter called as child grid). The cell size of the child grid $d_{cell}$ equals to $d_{CELL}/n_{sub}$, and $n_{sub}$ is an odd number bigger than 1 in this study to guarantee the symmetry of the simulated dendrite structure. That is, there are $n_{sub}^{3}$ child-grid cells in a parent grid cell. Similarly, there are several variables associated with the child-grid cell, they are $f_{s}$, the solid fraction, and $s$, the status of the cell ($s=0$: solid, and $s=1$: liquid).

In this model, the decentered octahedron growth algorithm (Fig. 1b) is implemented on the child grid to reproduce the movement of the solid-liquid interface. As Fig. 1a shows, the octahedral envelopes are located at the interfacial cells and constitute the solid-liquid interface (see the white dashed line). The six half-diagonals of these octahedral envelopes correspond to the crystallographic orientation of the dendrite structure, which implicitly incorporates the anisotropic surface energy. More details of this algorithm can be referred to Rappaz [14] and will not be repeated here. During the simulation, the child grid will be updated with the advancement of the solid-liquid interface: (1) the child grid in the parent-grid cell, which has fully solidified and has no envelopes, is removed; (2) the potential interfacial cell turns into an interfacial cell once at least one of its child-grid cells is captured; (3) new child grid will be generated at the liquid cells that are adjacent to the new interfacial cell.

The “handshaking” between the parent grid and the child grid is illustrated in Fig. 1c, where algorithms above the red dashed line are exerted on the parent grid, and those below the line are applied on the child grid. During the simulation, (1) the octahedral envelopes grow up by the movement calculated according to the interfacial solute concentration ($C_{l}$) and interface curvature ($\kappa$) at the parent grid scale; (2) envelopes capture their neighboring cells, updating the solid fraction ($f_{s}$) in the child grid; (3) the solid fraction ($F_{s}$) and the solute concentration ($C_{l}$) in the parent grid are updated; (4) the diffusion is calculated at the parent grid. These four procedures link the parent grid with the child grid and would be repeated till the end of the simulation. The detailed algorithms are introduced in the following sections.

In this model, the driving force for the solid-liquid interface movement is the deviation of the local liquid composition ($C_{l}$) from the equilibrium composition ($C_{l}^{eq}$). Within one time step ($\Delta t$), the change of solid fraction $\Delta F_{s}$ in the parent-grid cell can be calculated as:

where $k$ is the solute partition coefficient, $C_{0}$ is the initial solute concentration, $T_{0}^{eq}$ is the equilibrium temperature for the initial composition, $T$ is the local temperature, $\Gamma$ is the Gibbs-Thomson coefficient, $m_{l}$ is the liquidus slope, and $\kappa$ is the interface curvature and calculated using the counting-cell method [27] in the parent grid. The 3D version of this method is used in the present model as Eq. 3.

where 81 is the number of cells with the distance from the interfacial cell being within $\sqrt{6}d_{CELL}$, $\sum F_{s}$ is the summed solid fraction of cells that belong to the same dendrite as the interfacial cell do among the 81 neighboring cells, 51 is the summed solid fraction when the interface is flat, i.e. $\kappa=0$, and $\alpha$ is a shape factor equal to 4.3899, which is derived by regressing the calculated radii ($1/{\kappa}$) on the predefined sphere radii.

Thus, the interface movement ($\Delta s$) is:

In this model, the tips’ movement of all envelopes in the parent-grid cell is defined to be the same as the interface movement calculated on the parent grid, i.e., $\Delta s$. In the end of this time step, the envelope size $L$ which is the distance from the envelope center to its surface becomes:

where $L_{0}$ is the envelope size in the beginning of this time step.

To simulate the evolution of the solid fraction in child-grid cells, in the present study, a solid fraction calculation method is designed and applied on the child grid. $d$ is the distance from the envelope center to the center of cell $Q$ along the normal direction of the envelope surface closest to cell $Q$, $L_{c}$ is the maximum size of new envelopes and equals to $\sqrt{2}d_{cell}$ according to the decentered octahedron growth algorithm. The solid fraction in cell $Q$ ($f_{s}^{Q}$) updates as follows:

• 1)

  $L<L_{c}$: as shown in Fig. 2a, the envelope is too small to engulf cell $Q$, so $f_{s}^{Q}=0$.

• 2)

  $L_{c}<L<d$: as shown in Fig. 2b, cell $Q$ is outside of the envelope but partially engulfed by the envelope, the solid fraction will be:

  $$f_{s}^{Q}=\frac{L-L_{c}}{d-L_{c}}$$

  (6)

• 3)

  $L\geq d$: as shown in Fig. 2c, cell $Q$ is captured by the envelope, then $f_{s}^{Q}=1$. At the same time, a new envelope hosted by cell $Q$ is generated according to the decentered octahedron growth algorithm.

If a cell is the neighboring cell of several envelopes, its solid fraction will be the maximum $f_{s}$ calculated as Eq. 6. After the completing the growth and capture process of all envelopes, the solid fraction in the parent-grid cells ($F_{s}$) can be obtained by averaging the solid fraction of its child-grid cells as Eq. 7.

As the dendrite grows up, the solute partition occurs at the interfacial cells. The solute concentration of the newly solidified child-grid cells will be $k\cdot C_{l}$. The rejected ($(1-k)C_{l}$ when $k<1$) or absorbed ($(k-1)C_{l}$ when $k>1$) is evenly added to or subtracted from the liquid part of this parent-grid cell and its neighboring parent-grid cells.

As shown in Fig. 1c, the solute diffusion model is defined on the parent grid. Usually, the diffusion coefficient in liquid is several orders of magnitude higher than that in solid, thus only diffusion in liquid is accounted here. In this study, the solute diffusion is modeled using the Fick’s law:

where $D$ is the effective diffusion coefficient and dependent on the cell status. In this model, the effective diffusion coefficient ($D_{AB}$) between cell $A$ and one of its Von-Neumann neighboring cell $B$ is defined as:

where $D_{l}^{A}$ and $D_{l}^{B}$ are the diffusion coefficients in cell $A$ and cell $B$, respectively, $F_{s}^{A}$ and $F_{s}^{B}$ are the solid fractions in cell $A$ and cell $B$, respectively. Based on this equation, the diffusion coefficient is 0 in the fully solidified region, $D_{l}$ in the liquid area, and ranges from 0 to $D_{l}$ in the interfacial region.

The time step size ($\Delta t$) is predefined before the simulation starts and equals to $(d_{CELL})^{2}/(4D_{l})$ in this study in order to guarantee the stability of the solute diffusion model. Within one time step, the calculation scheme is as follows:

• 1)

  Loop over all interfacial parent-grid cells and calculate the interface movements.

• 2)

  Loop over all envelopes to complete the growth and capture processes and calculate the solid fraction in child-grid cells, then update the solid fraction in parent-grid cells.

• 3)

  Update the child grid as presented in Section 2.1.

• 4)

  Loop over all parent-grid cells with solid fraction changed in this time step to update its solute concentration.

• 5)

  Calculate the diffusion and then update the temperature field if necessary.

The dendrite growth in Fe-0.6 wt.%C alloy melt is simulated in this section. The critical properties of the material are listed on Table. 1.

Yasuda et al. [4] observed the dendrite growth process in Fe-0.58 wt.%C alloy melt using the x-ray imaging technique and presented the process in the form of video, which is used to directly validate our model. In the video, from $t=-37$ s to $t=-9$ s and before the $\delta\rightarrow\gamma$ phase transformation occurs, the dendrite (red ellipse in Fig. 8, and $\varphi=64^{\circ}$ ) grows steadily from 0.20 mm to 1.15 mm. In this study, the evolution of the dendrite trunk length is extracted and plotted in Fig. 10, from which we can see that growth velocity ($V_{m}$) is relatively constant around 0.0337 mm/s. As described in the reference [4], during the dendrite growth process, the cooling rate ($R_{c}$) is 0.17 K/s and there is a temperature gradient ($G_{T}$) along the vertical direction. However, the exact temperature gradient is not given in their article. Nevertheless, since the growth velocity ($V_{m}$) is relatively constant, it is reasonable to consider the undercooling at the dendrite tip to be invariable, which then leads to the equation $V_{m}\cdot sin(\varphi)\cdot G_{T}=R_{c}$. Therefore, the temperature gradient is approximated to be $5.6\times 10^{3}$ K/m.

The simulation is conducted on a $80\times 250\times 300$ mesh with the parent-grid cell size $d_{CELL}=4$ $\rm\mu m$ and the child-grid cell size $d_{cell}=d_{CELL}/3$. The initial temperature at the bottom of the simulation domain is 1762.75 K, i.e. the undercooling is 2 K, which is determined using the measured growth velocity ($V_{m}$) based on the K-G-T model. Additionally, the cooling rate and the temperature gradient are the same as those in the experiment, and the time step size is calculated to be 2 ms. At the beginning, a dendrite seed with the same crystallographic orientation as that of the dendrite highlighted in Fig. 8b is placed on cell (40, 150, 0). During the simulation, the dendrite is not allowed to cross the domain boundary in this model. This simulation is carried out on a desktop equipped with Intel core i7-9700 and costs about 24 hours.

The evolution process of the dendrite is presented in Fig. 9. The dendrite exhibits a branchless shape at the very beginning and then develops its secondary arms (Fig. 9a, b, e and f). Afterwards, the growth of the dendrite trunk at the left side is blocked by the domain boundary. The dendrite trunk stop growing longer. Some secondary arms on it that are parallel to the dendritic trunk pointed by the red arrow develop the tertiary arms (see Fig. 9c,d, g and h), resulting in the shape of a new dendrite trunk. In Yasuda’s video [4], although there is no boundary restricting the dendrites’ growth, the left branch of the dendrite (highlighted by red circle) is blocked by the neighboring dendrite. In this branch, one of the secondary dendrite arms develop the tertiary arms, and finally become a new dendrite trunk (see Fig. 8b), which is similar to the phenomenon in our simulation. Besides, as the cross-section views in Fig. 9e, f, g and h show, when the dendritic trunk grows towards the top of the domain, dendritic arms at the bottom coarsen continuously. These phenomena are also observed in Yasuda’s video [4]. Furthermore, the length of the dendritic trunk (pointed by a red arrow in Fig. 9) at different times are also extracted and plotted in Fig. 10. Obviously, the simulation result is in good agreement with the experimental data. According to the growth kinetics represented by Eqs. 1, 2 and 3, the accurate evaluation of the growth depends on the evaluation of the interface curvature, which needs the accurate capture of the solid-liquid interface. Therefore, the good agreement in the trunk growth curve also indicates the good capture of the solid-liquid interface in the present model.

In additive manufacturing, especially the Electron Beam Selective Melting (EBSM) of Ni-based superalloy, the alloy solidifies with a significant segregation in the inter-dendritic regions despite of the extremely high temperature gradient ($\approx 10^{5}-10^{7}$ K/m) and cooling rate ($\approx 10^{3}-10^{5}$ K/s). This process usually comes with the formation of precipitates and the incubation and propagation of hot cracks, which would affect the mechanical performance of the as-built products. Simulating the dendrite growth process during the solidification process would help us understand the physical mechanism, and then provide guidance on how to tailor the microstructure evolution during the AM process.

Single track is the fundamental unit of additive manufacturing. In this study, to demonstrate the model’s validity in AM condition, a single-track scanning experiment on an Inconel 718 substrate is carried out in a self-developed EBSM machine. The substrate was preheated at the beginning with the defocused electron beam. The scanning with focused electron beam started when the temperature of the substrate reached 1123 K. The beam power is 400 W, and the scanning speed is 0.5 m/s. The sample was cut along the plane perpendicular to the scanning direction in order to show the cross section of the track. The cross section was polished and then etched for microstructure observation in a scanning electron microscope (SEM). Meanwhile, energy dispersive X-ray spectroscopy (EDS) analysis was conducted to characterize the chemical concentration field. The single track has a width of 702 $\rm\mu m$ and a depth of 219 $\rm\mu m$. In Fig. 11e, three regions that are located near the center plane are chosen for the comparison with the simulation results using the presented dendrite growth model. The distance between the regions and the bottom of the melting pool are 12 $\rm\mu m$, 48 $\rm\mu m$ and 84 $\rm\mu m$, respectively. In the melting pool, the dendrites prefer to grow towards the center of the melting pool, making the dendrite structure at the upper part of the melting pool extremely complex (see Fig. 11e). Therefore, the selected regions are all at the lower part of the melting pool, where the critical indicator, i.e., the primary dendrite arm spacing (PDAS), is measurable. The PDASs in these regions are measured and plotted in Fig. 13a.

A thermal-fluid model developed by Yan et al. [34] is used to reproduce the melting and solidification process. In this model, the major factors including the fluid flow, Marangoni effect, metal evaporation, and recoil pressure are incorporated, so that a reliable thermal field can be obtained, which has been validated in our previous work [35, 36]. In this simulation, the initial temperature of the substrate was set as 1123 K instead of modelling the preheating process. The other scanning parameters are exactly the same as those used in the experiment. The beam diameter was set as 700 $\rm\mu m$ by calibrating against the experimental result. The material properties are listed in Table 2. The simulation result is shown in Fig. 11. Obviously, the comparison of the cross sections perpendicular to the scanning direction in Fig.11b indicates that the scanning process is well reproduced by the thermal-fluid model. According to the distances between the selected regions and the bottom of the melting pool (Fig. 11e), the corresponding solidification conditions (temperature gradient and cooling rate) are extracted from the thermal-fluid simulation (Fig. 11c and d). The detailed data are shown in Table 3. These parameters are then used in simulating the dendrite growth in the three regions, respectively.

In this study, the dendrite growths in the three regions are carried out using the temperature gradients and cooling rates listed in Table 3. Note that the temperature gradient and cooling rate are constant in the dendrite growth model. Since Nb is one of the most important segregation elements in Inconel 718, and the Nb-rich precipitates in the inter-dendritic regions significantly affect the mechanical performance of the as-built products, the material is treated as Ni-Nb binary alloy following [22, 37, 38]. The concentration of Nb is 5 wt.%. The detailed physical properties are listed in Table. 4. The simulation domain is a $60\times 60\times 100$ mesh with $d_{CELL}=0.2$ $\rm\mu m$ and $n_{sub}=3$. Before the simulations start, a number of dendrite seeds (about 420) with the same crystalline orientation are randomly placed on the bottom of the simulation domain. The distance between these dendrite seeds is restricted below 0.6 $\rm\mu m$. The initial temperature at the bottom of the simulation domain is set as the liquidus temperature of Inconel 718.

As shown in Table 3, the solidification velocities are high so that the solute trapping cannot be ignored. In this study, the solidification velocity ($V$)-dependent solidification parameters (partition coefficient $k_{v}$ and the liquidus slope $m_{v}$) are calculated according to Eq. 16 (Aziz’s model [39]) and Eq. 17 [30], and adopted in the simulations. In Eq. 16, $V_{D}$ is the interface diffusion velocity, which is material-dependent. Up to date, there is no precise measurement of $V_{D}$ for Inconel 718. Ghosh et al. [40] fitted the simulated solidification velocity-partition coefficient relation with Aziz’s model, thereby determining $V_{D}=0.62$ m/s. This result is consistent with the deduction ($V_{D}<1$ m/s) made by Tao et al. [3] based on their experimental results. Therefore, in the study, $V_{D}=0.62$ m/s, and then the corresponding $k_{v}$ and $m_{v}$ for the three regions are calculated (see Table 3) and adopted in the dendrite growth simulations.

Fig. 12 shows the simulated dendrite structures when the tip undercooling is under the steady state. Obviously, the number of dendrites decreases from region 1 to region 3. It is critical to note that the number of dendrites is much smaller than that of the dendrite seeds placed on the domain bottom, which is due to the growth competition among dendrites. In this study, the simulated PDAS are calculated by $\lambda_{1}=\sqrt{A/N_{d}}$ ($A$ is the area of the XY cross section, and $N_{d}$ is the number of remaining dendrites observed from the domain top). As shown in Fig. 13a, the simulated PDASs are quite close to the measured results. The tip undercoolings obtained when the dendrite growth reaches the steady state are plotted in Fig. 13b. The theoretical values (considering the $V$-dependent $k_{v}$) calculated by [40] are also plotted for comparison. Both the theoretical values and the simulated values indicate that the higher solidification velocity comes with a higher tip undercooling. Although the simulated results are all lower than the theoretical predictions, the differences between the simulation results and the theoretical predictions are not that big ($V=0.02$ m/s: 6.5 K or 27.7%; $V=0.05$ m/s: 3.9 K or 13.7% and $V=0.075$ m/s: 2.1 K or 7%). That is, the simulation results are fairly close to the theoretical predictions. The differences are mainly the impact of the discretization in the model. Especially, since $d_{CELL}=0.2$ $\rm\mu m$, the temperature change $\delta T$ between neighboring cells along the $Z$ direction are 1 K ($V=0.02$ m/s), 0.2 K ($V=0.05$ m/s) and 0.08 K ($V=0.075$ m/s), respectively. Obviously, the smallest $\delta T$ in region 3 results in the smallest difference as shown in Fig. 13b. Additionally, in comparison with the phase field model in [40], the present model shows a much better consistency with the theoretical prediction when the solidification velocity is higher than 0.03 m/s.

As Nb-rich precipitates significantly affect the mechanical performance of the as-built parts, the simulated Nb distribution is compared with the measured result. Firstly, both the simulation result and experimental observation show that Nb prefers to concentrate at places where the dendrites compete for growth space fiercely. For example, the regions highlighted by red circles in Fig. 14a, and the middle region in Fig. 14c. At these regions, the Nb concentrations are significantly higher than those at ordinary inter-dendritic regions. Secondly, the Nb concentrations at representative points (see Fig. 14b) are listed in Table 5. Point 3 and 5 are located at the dendrite trunks, the Nb concentrations are similar to the simulation results at the dendrite trunks. Thirdly, in the simulation result and the experimental result, there are some small regions with very high Nb concentration, such point 1, 2 and 4 in Fig. 14b, and red regions in Fig. 14c. These regions are all located at the fierce-competition places. Obviously, the measured Nb concentration in point 1, 2 and 4 are all higher than the maximum value in the simulation result, which is due to the further element concentration caused by the formation of Nb-rich precipitates (see the corresponding points in Fig. 14a). Detailed investigation into this issue is not in the scope of the present work.