GPU-Based Cellular Automata Model for Multi-Orient Dendrite Growth and the Application on Binary Alloy

Abstract

To simulate dendrite growth with different orientations more efficiently, a high-performance cellular automata (CA) model based on heterogenous central processing unit (CPU)+ graphics processing unit (GPU) architecture has been proposed in this paper. Firstly, the decentered square algorithm (DCSA) is used to simulate the morphology of dendrite with different orientations. Secondly, parallel algorithms are proposed to take full advantage of many cores by maximizing computational parallelism. Thirdly, in order to further improve the calculation efficiency, the task scheduling scheme using multi-stream is designed to solve the waiting problem among independent tasks, improving task parallelism. Then, the present model was validated by comparing its steady dendrite tip velocity with the Lipton–Glicksman–Kurz (LGK) analytical model, which shows great agreement. Finally, it is applied to simulate the dendrite growth of the binary alloy, which proves that the present model can not only simulate the clear dendrite morphology with different orientations and secondary arms, but also show a good agreement with the in situ experiment. In addition, compared with the traditional CPU model, the speedup of this model is up to 158×, which provides a great acceleration.

1. Introduction

Dendrite is an important component of solidification structure, whose morphological characteristics directly affect the performance of an alloy. Therefore, studying the growth of dendrite is helpful for researchers to optimize the process parameters and improve the performance of alloys. Experimental methods for obtaining the microstructure during solidification have been proposed, such as the electronic probe [1], synchrotron X-ray [2,3], etc. that have been developed, but cannot be applied to product at a large scale due to the expensive cost and harsh environment. Moreover, these methods are unable to show the dynamic evolution process of dendrite growth.

Over the past few decades, a considerable number of numerical models have grown up around the theme of simulating the evolution process of dendrite morphology during solidification, such as level set [4], monte carlo (MC) [5], cellular automata (CA) [6,7,8], phase field (PF) [9], etc. Among them, the phase field model (PF) and CA model are the most widely used. Due to complex physical equations and massive calculations, the PF model is only applied to small-scale simulations. Whereas the CA model can perform large-scale simulations and reveal the evolution of micro-mesoscopic dendrites, such as the transformation of columnar crystals into equiaxed crystals (CET) [10], the behavior of dendrite deflection in fluids [6,7,8], and the formation of silicon facet dendrite [11]. Influenced by grid layout and capture rules, the traditional CA model can only simulate dendrites that grow along the x axis or 45° with the x axis. However, during the actual solidification process of the alloy, dendrites grow with different preferred orientations, which is caused by the solute profile and can lead to macro-segregation in the cast product [12]. Therefore, it is very important for the model to be able to simulate the dendrite with different orientations. Researchers have proposed various methods to achieve multi-orient dendrite simulation, such as defining block cells [13], establishing random ZigZag neighbor cell capture rules [14], etc. These methods can only simulate a dendrite with the same orientation at the same time due to complex algorithms. Hence, Rappaz and Gandin developed the decentered square algorithm (DCSA) [15], which can simulate dendrites with different orientations simultaneously, and it has been widely used to simulate the growth of columnar and equiaxed dendrites [16,17,18,19]. However, these models suffered from velocity anisotropy, which means velocities along different orientations are asymmetrical and have multiple interfaces. Then, Luo et al. [20] modified the velocity of the interface cell through the local solute equilibrium and improved the phenomenon of multiple interfaces. Wang et al. [21]. determined the length of the square half diagonal by the preferred orientation to ensure sharpness of the interface and introduced GF to reduce velocity anisotropy.

However, the DCSA rule dynamically determines the capture position and condition by tracking the dendrite tip, which is very complex. Massive calculations of the CA model incorporating the DCSA rule make large-scale simulations very time-consuming. It is, therefore, necessary to speed up the calculation. Some methods have attempted to improve calculation efficiency, such as adaptive meshes [22,23], parallelizing the computer program using MPI technology [24], and parallel computing technology based on serial arithmetic [25]. These accelerators are CPU-based, so the speedup is not good enough for the number of cores, being limited by the hardware architecture. Unlike the CPU, GPU is characterized by integrating thousands of computer cores inside, which makes it the most effective way to achieve acceleration. With the development of GPU in general purpose computing, it has been used in various fields [26,27]. Over the past two decades, the GPU-based acceleration method has become a hotspot in computational and material science [28,29,30] and has made great progress in dendrite growth simulation [31,32,33,34]. However, most of these contributions focused on the PF model, which is still unable to simulate dendrite growth at large scale, even when accelerated by using GPU.

This paper proposes a high-performance CA model incorporating the DCSA capture rule to simulate dendrite growth with different orientations more efficiently. To make full use of the hardware resources and improve the calculation efficiency, parallel algorithms and a task scheduling scheme using multi-stream are proposed. Compared to the traditional CPU-CA model, this model can achieve great acceleration with a speedup of 158×. In addition, the steady-state tip velocity predicted by this work shows great agreement with the analytical value of the LGK model. The simulation results of the binary alloy by this model not only show dendrites with different orientations, but also agree well with the experimental result in situ. This model will be a promising tool for studying the dendrite growth during the large-scale solidification.

The following chapters will cover the establishment, solution, and verification of the model in detail.

2. Materials and Methods

2.1. Description of Numerical Model

2.1.1. Heat Transfer Model

Equation (1) is used to describe heat transfer during alloy solidification.

$$p{c}_p\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left(\lambda \frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left(\lambda \frac{\partial T}{\partial y}\right)$$

(1)

where $T$ is the temperature, $f_s$ is the solid fraction, $\rho$ is the density, $c_p$ is the equivalent specific heat, and $L$ is the latent heat.

2.1.2. Solute Distribution Model

The solute distribution model consists of solute redistribution at the interface and solute diffusion across the entire domain. For one thing, the solute redistribution is determined by Equations (2) and (3) according to local equilibrium.

$$C_s=k{C}_l$$

(2)

$$dC=C_l\left(1-k\right)\mathsf{\Delta }{f}_s$$

(3)

where $dC$ is the residual solute expelled by the interface cell, and $k$ is the coefficient of equilibrium solute partition ($k<1$ in the present work).

First, the expelled solute will be discharged into the remaining liquid of the current cell. Then, the current $C_l$ will be compared to the equilibrium solute concentration $C_l{}^{*}$; the part beyond $C_l{}^{*}$ will be discharged to adjacent liquid cells. $C_l{}^{*}$ is calculated by Equation (4).

$$C_l^{*}=C_0+\frac{T-T_0}{m_l}+\frac{\Gamma Ka\left(\widehat{n}\right)}{m_l}$$

(4)

where $T$ is the actual temperature, $T_0$ is the initial temperature, $C_0$ is the initial solute concentration, $m_l$ is the liquidus slope, and $a\left(\widehat{n}\right)$ is the function of interface anisotropy.

If the current cell is completely solidified, then all the excluded solute will be discharged to the adjacent liquid cells. How much solute will be discharged into an adjacent liquid cell is designed by weight, as shown in Equation (5).

$$w_i=\frac{C_i^0-C_l}{{{\displaystyle \sum }}_{j=0}^N(C_j^0-C_l)}$$

(5)

where $C_i{}^0-C_l$ is the solute difference between the interface cell and the adjacent liquid cell $i$, $N$ is the number of liquid neighbours.

For another, the solute diffusion coefficient in liquid is three to four orders of magnitude larger than that in solid, so diffusion in solid phase is neglected in this paper. Equation (6) shows the diffusion of a solute in liquid and interface.

$$\frac{\partial C_l}{\partial t}=\frac{\partial }{\partial x}\left(D_i\frac{\partial C_l}{\partial x}\right)+\frac{\partial }{\partial y}\left(D_i\frac{\partial C_l}{\partial y}\right)$$

(6)

where $D_i$ is the solute diffusion coefficient, which is governed by Equation (7).

$$D_i=f_s{D}_s+\left(1-f_s\right)D_l$$

(7)

where $D_l$ and $D_s$ are the diffusion coefficient in liquid and solid, respectively.

2.1.3. CA Model

In the traditional CA model, only fully solidified cells ($f_s=1$) can capture their liquid neighbors. This capture rule is severely affected by the grid layout, so it can only simulate dendrites growing along the x axis and 45° with the x axis. However, the growth orientation of dendrite is chaotic in the actual solidification process of the alloy. In order to simulate dendrite growth closer to reality, the decentered square rule for capturing has been adopted by the present CA model. As shown in Figure 1, the key to the improved algorithms is dynamically calculating the nucleation position and nucleation conditions, according to the position where the interface cell (parent cell) reaches the liquid cell. All cells in the entire computing domain are initialized to liquid ($f_s=0$). As the temperature decreases, the one who receives the nucleation conditions will receive the chance to nucleate. The nucleation probability is calculated by Equation (8) combined with random probability.

$$\frac{dn}{d\mathsf{\Delta }T}=\frac{n_{max}}{\sqrt{2}\exp \left[-\frac{1}{2}\left(\frac{\mathsf{\Delta }T-{\mathsf{\Delta }T}_N}{{\mathsf{\Delta }T}_{\sigma }}\right)\right]}$$

(8)

where $\mathsf{\Delta }T$ is the undercooling, $n_{max}$ is the maximum density of nuclei given by the integral of undercooling (from 0 to $\mathsf{\Delta }T$) of Equation (8), $\mathsf{\Delta }{T}_N$ is the average nucleation undercooling, and $\mathsf{\Delta }{T}_{\sigma }$ is the standard deviation. Once a cell is nucleated, its solid fraction will be changed from 0 to 1, and a decentered square with an angle $\theta$ between diagonal and x axis will be located at its center. Additionally, four corners of this square penetrate the nearest liquid neighbors, capturing them in the corners and changing them into interface cells. Captured cells will inherit the growth properties of the parent cell and grow along the diagonal to capture its liquid neighbors until it is completely solidified.

In order to simulate the dendrite morphology in detail, the local level rule method is used to calculate the increment of the solid fraction [35], wherein $\mathsf{\Delta }{f}_s$ of each time step is determined by Equation (9). The length of the primary arm of the dendrite is updated by Equation (10).

$${\mathsf{\Delta }f}_s=\frac{\left(C_l^{*}-C_l\right)}{\left(C_l^{*}\left(1-k\right)\right)}GF$$

(9)

$$l\left(t+\mathsf{\Delta }t\right)=l\left(t\right)+{\mathsf{\Delta }f}_s{l}_{max}$$

(10)

$l_{max}$ is the maximum half-diagonal length determined by Equation (11).

$$l_{max}=l\left(\frac{\mathsf{\Delta }x}{\max \left(\sin \theta ,\cos \theta \right)}\right)$$

(11)

$\mathsf{\Delta }x$ is the size of the cell, and $GF$ [36] is the geometry factor related to the state of neighbors, which is introduced to avoid multiple interfaces and is limited to no more than 1.

$$GF=\min \left\{1,G({{\displaystyle \sum }}_{m=1}^4{S}_m^{\mathrm{I}}+\sqrt{2}{{\displaystyle \sum }}_{m=1}^4{S}_m^{\mathrm{II}})\right\}$$

(12)

where $G$ is an adjustable factor, $s_m{}^{\mathsf{{\rm I}}}$ and $s_m{}^{\mathsf{{\rm I}}\mathsf{{\rm I}}}$ are the state of the nearest neighbor and the second nearest neighbor, respectively, and both of them have two values, as shown in Equation (13).

$$S_m^{\mathrm{I}},S_m^{\mathrm{II}}=\{\begin{array}{c}0, f_s<1\\ 1, f_s=1\end{array}$$

(13)

2.2. Parallel Solver Based on GPUq

To simulate the morphology of dendrite with clequeqear secondary arms, the cell size is generally divided into 1~2 $\mathsf{\mu }\mathrm{m}$. There will be huge calculations when conducting large-scale simulations, which is time-consuming using CPU. In addition, the purpose of the simulation is to optimize the process parameters by performing many groups of experiments. Therefore, it is necessary to accelerate the calculation. This paper designs a GPU-based parallel solver for acceleration.

2.2.1. GPU-CA Framework

The GPU-CA frame used in this work is shown in Figure 2, in which each slice has the cell state of the cross-section corresponding to this moment.

In this heterogeneous framework, the CPU is responsible for allocating memory, initializing, and transferring the initial state to the global memory on the GPU via the PCIe bus. Then, all the threads on the GPU execute instructions simultaneously in single-instruction-multiple-data (SIMD) mode until the required calculations are completed. Finally, data are transferred from GPU memory to CPU memory via the PCIe bus for post-processing. This work uses Pascal GP100 GPU (embedded with 3584 CUDA cores) as the accelerator, where millions of threads can execute instructions simultaneously, and each thread corresponds to a cell.

2.2.2. Implementation by CUDA C

This work adopts CUDA as the programming model, because it allows us to execute applications on heterogeneous computing systems by simply annotating code with a small set of extensions to the C programming language. In addition, the CUDA programming model provides the exposed memory hierarchy, in which global memory is the most widely used, because it has a large memory space. This work chose global memory and, based on it, designed the structure of array (SOA) to organize data, because the SOA is more liable to maximize the efficiency of accessing the global memory.

• Task scheduling scheme

CUDA provides streams to achieve concurrency between kernels, where kernels are executed sequentially in the same stream, simultaneously in different streams. To increase parallelism, this paper divides the problem into eight small tasks, each corresponding to a kernel function, as shown in

Table 1. Task assignment in GPU-CA model.

Kernel Functions	Computation Task
$\mathrm{capture}$	capturing neighbors
$\mathrm{calD}$	solute diffusion coefficient
$\mathrm{solute}\_\mathrm{Dif}$	solute diffusion
$\mathrm{schange}$	solute redistribution in interface
$\mathrm{get}\_\mathrm{T}$	temperature distribution
$\mathrm{calB}$	equilibrium solute
$\mathrm{growth}$	velocity and the arm length
$\mathrm{backup}$	storing data for the next slice

.

All kernels are related, but at some point, some kernels are independent. For example, before kernel $\mathrm{solute}\_\mathrm{Dif}$ is scheduled, kernels $\mathrm{capture}$, $\mathrm{schange}$, $\mathrm{calD}$ are independent. To solve the waiting problem among independent tasks, this paper puts all kernels into two streams based on the kernel dependency, as shown in Figure 3. It is noted that all operations in the non-default stream are non-blocking with respect to the host thread. Thus, we need to synchronize the host with operations running in a stream. The overall flow of the numerical simulation is shown in Figure 3, where kernel function $\mathrm{solute}\_\mathrm{Dif}$ depends on values of $C_L$ and $D$ calculated by $\mathrm{calD}$ and $\mathrm{schange},$ respectively. In this case, the host has to wait until the $\mathrm{calD}$ and $\mathrm{schange}$ finish their calculation before $\mathrm{solute}\_\mathrm{Dif}$ is scheduled. This work uses the $\mathrm{cudaDeviceSynchronize}\left(\mathrm{void}\right)$ function to achieve synchronization between the host and operations running in a stream.

2

Parallel algorithms

Global memory can be accessed on the device from any SM throughout the application’s lifetime. When multiple threads write data to the same address at the same time, there will be a ‘data race’ that can lead to an undefined error. The CA model in this work adopted the Moore neighbor, in which each cell has eight neighbors. Therefore, there may be more than one neighbor discharging solute into a liquid cell or capturing it at the same time. Mapping to the programming model, there are multiple threads that attempt to modify the data in the same address, introducing a ‘data race’. This paper avoided this phenomenon by preventing the cell from ‘actively’ discharging solute to neighboring cells or capturing it. For example, the solute redistribution algorithm adopts a kernel function and a device function to turn ‘active’ to ‘passive’, as shown in Figure 4. There are two steps to complete the solute redistribution process. The pseudo-code is listed in Algorithm 1.

First, calculate how much solute will be discharged to each liquid neighbor and store it in the intermediate variables, lines 2–7 in Algorithm 1.

Then, read the data stored by the first step from the neighbors and sum it to the solute of the current cell changed, lines 8–11 in Algorithm 1.

Algorithm 1 parallel solute redistribution algorithm

1 $id: \mathrm{assign}\mathrm{thread}\mathrm{to}\mathrm{cell}$

2 $\mathrm{for}\mathrm{all}\mathrm{interface}\mathrm{cells}\mathrm{do}$

3  $dC\leftarrow \mathrm{diacharged}\mathrm{solute}$

4  $\mathrm{if}(dC>0)$

5   $d{C}_i\leftarrow \mathrm{store}\mathrm{the}\mathrm{solute}\mathrm{dischared}\mathrm{to}\mathrm{cells}$

6  $\mathrm{end}$

7 $\mathrm{end}$

8 $\mathrm{for}\mathrm{all}\mathrm{liquid}\mathrm{and}\mathrm{interface}\mathrm{cells}\mathrm{do}$

9  ${\mathrm{temp}\leftarrow \mathrm{read}\mathit{dC}}_i\mathrm{from}\mathrm{neighbors}$

10  $C_l\leftarrow C_l+\mathrm{temp}//\mathrm{add}\mathrm{temp}\mathrm{to} C_l$

11 $\mathrm{end}$

This algorithm can avoid the ‘data race’ and make all cells in the calculation domain independent.

3. Results and Discussion

3.1. Validated by the LGK Model

This model was validated by comparing the steady velocity of a free dendrite crystal in an undercooled melt with the analytical value of the LGK model. This work defines the velocity of the dendrite tip through the cell size and time interval as the cell stays in the interface state, and the steady state of dendrite growth is determined as the solute at the boundary opposite the dendrite tip reaches 1.01 times the initial value [37,38]. The physical parameters used in this paper are shown in

Table 2. Physical properties [39,40].

Property and Symbol	Fe–0.6C	Fe–5.3Si
Initial composition,$C_0\left(\mathrm{wt}.\%\right)$	0.6	5.3
Liquidus temperature $T_l\left(\mathrm{K}\right)$	1763.37	1732.87
Liquidus slope, $m_l\left({\mathrm{K}\%}^{-1}\right)$	−80	−7.6
Solute partition coefficient, $k_0$	0.34	0.77
Solute diffusion coefficient in liquid $D_l\left({\mathrm{m}}^2·{\mathrm{s}}^{-1}\right)$	2 × 10−9	8.0 × 10−4exp (−29,943.23/T)
Solute diffusion coefficient in solid $D_s\left({\mathrm{m}}^2·{\mathrm{s}}^{-1}\right)$	5 × 10−10	8.0 × 10−8exp (−29,943.23/T)
Gibbs-Thomson coefficient,$\left(\mathrm{m}·\mathrm{K}\right)$	1.9 × 10−7	1.9 × 10−7

.

This test was performed in a square domain with a size of 400 μm × 400 μm, which is uniformly divided into 400 × 400 cells with the size of 1 μm × 1 μm. At the beginning of solidification, a nucleus with a preferred orientation of 0° was placed in the center of the domain.

Figure 5 shows the evolution of the transient velocity of the dendrite tip as a function of time calculated by the present CA model under the undercooling of 8K and the steady velocity calculated by the LGK model. It can be seen that tip velocity of the dendrite decreases dramatically in the transient growth period; then, the decreasing tendency becomes mild in the steady growth period. Finally, the velocity is stable at 56.33 μm/s, which is very close to the analytical value of 55.90 μm/s calculated by the LGK model.

3.2. Model Capability

This model is applied to simulate the dendrite growth of binary alloy. Firstly, a group of dendrites of Fe–0.6C alloy are simulated to evaluate the ability to simulate the dendrites with different orientations. Secondly, it is applied to simulate multi-orient dendrites’ growth of Fe–5.3Si alloy [41], and the simulated result is compared with the in situ observation. Finally, the performance of the parallel solution is analyzed by comparing the calculation time with the traditional CPU calculation.

3.2.1. Single Dendrite of Fe–0.6C Alloy

A group of dendrites of Fe–0.6C with different orientations are simulated in this part, where the orientation ranges from 0° to 90° due to the four-fold symmetry of the dendrite. Simulations are also performed on a square domain with a size of 400 μm × 400 μm that is divided into $400\times 400$ cells. Figure 6 shows the simulation results at 0.6 s with a constant undercooling of 8K, which shows that the present model can simulate clear dendrite morphology with not only different orientations but also second dendrite arms.

In addition, as shown in Figure 7, the primary arm of dendrites with different orientations of 0°, 30°, and 60° are almost the same length, indicating that growth velocities with different orientations are in good symmetry.

3.2.2. Multi-Dendrites of Fe–5.3Si Alloy

The development of synchrotron X-ray technology has enabled the observation of solidification in metallic alloys, which provides a powerful tool to verify the model. Yasuda et al. [41] performed in situ observation of Fe–5.3Si wt.% alloy, which shows a clear morphology of dendrite, including secondary dendrite arms. To validate this model, this paper simulated multi-orient dendrites’ growth of Fe–5.3Si alloy under the same experimental condition as the in situ observation, which set the nuclei positions and orientations corresponding to the in situ experimental observation before the simulation began. Figure 8a,b depict the simulated dendrite morphology and solute profile, respectively.

Simulation results show that the bottom dendrite is less developed than the upper because there are more dendrites at the bottom than at the upper, which agrees well with the in situ observation by synchrotron X-ray imaging [41]. For one thing, the growth of many dendrites will make the surrounding solute enriched, which will restrict the development of dendrites. For another, the more dendrites, the less space for each of them. In addition, Figure 8b shows the segregation between dendrites and within a dendrite, which satisfies well the non-equilibrium solidification theory [42]. Therefore, this model can simulate the growth of the dendrite in the actual casting process and reveal the phenomenon of segregation.

3.2.3. Acceleration Performance

The block size configured when launching a kernel function means a lot to performance because of its impact on latency hiding, memory efficiency, and occupancy, etc. A group of numerical experiments were carried out to find the optimal execution configuration according to the grid and block size guidelines [43]. Further, parts of the representative numerical simulation results are shown in

Table 3. Time elapsed with different kernel configuration.

Kernel Configuration	Time Elapsed (s)
(128,1)	37.94
(128,2)	37.64
(128,4)	38.85
(64,4)	37.59
(64,2)	37.53
(62,1)	38.13
(32,8)	37.70
(32,4)	37.23
(32,2)	37.85

, which shows that the configuration of (32,4) can obtain optimal performance. Accordingly, the following numerical experiments are performed under this configuration.

Speedup, the ratio between the CPU and GPU computational time, is used to evaluate the performance of the present GPU-CA model. The same numerical simulations were performed both on the Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40 GHz and Pascal GP 100 GPU. All simulations are carried out on the assumption that solidification has been started with 0.6 s, which equals 600,000 steps in the case of $\mathsf{\Delta }t=0.0001\mathrm{s}$. Figure 9 shows the computation time on the CPU and GPU as well as the speedup, which shows that this model can achieve great acceleration. In addition, the speedup increases with the number of grids, which is up to 158×, but the tendency to increase becomes low from a certain point. This can be explained by limited hardware resources, and before it is fully utilized, performance will be significantly improved when the grid number is increased. Once occupancy reaches its maximum, performance may be limited by additional scheduling overhead when increasing the grid number again.

4. Conclusions

This paper aims to develop a high-performance CA model incorporating the decentered square capture rule to simulate dendrite growth with different orientations more efficiently. The calculation efficiency of the CA model based on GPU has been improved by the proposed parallel algorithms adapted to the many-core architecture of GPU and the task scheduling scheme using multi-stream, and it is validated by comparing the calculated value with the analytical value by the LGK model. By applying this model to simulate the single dendrite of the Fe–0.6C alloy in different orientations and the multi-orient dendrite growth of the Fe–5.3Si alloy, the following conclusions can be drawn:

(1)

The steady dendrite tip velocity calculated by this CA model agrees well with the analytical LGK model.

(2)

The present model can simulate dendrite morphologies with a random orientation of the Fe–0.6C alloy and can maintain velocity symmetry under different orientations.

(3)

The simulation result of Fe–5.3Si not only matches well with the in situ experiment but can also reveal the segregation existing between and within dendrite. This model can be used to simulate multi-dendrite growth in actual casting.

(4)

Compared to traditional CPU calculation, this work can achieve noticeable acceleration, and the speedup increases with the number of grids, which is up to 158×.

This GPU-based model can greatly accelerate the calculation, which will be a promising tool in the field of studying dendrite growth during large-scale solidification.