A High-efficiency Virtual Submesh Cellular Automata Method for Solidification Simulation with Low Mesh Anisotropy

Abstract

The cellular automata (CA) method has been widely used in microstructure simulation during solidification, but it is time-consuming when used for multiple orientations and large-scale computation. In this work, a novel high-efficiency virtual submesh cellular automata (VISCA) method is proposed to reduce the mesh-induced artificial anisotropy. This method conceptually divides the mesh into a series of submeshes by giving a unique index to each CA cell. The proposed VISCA method results in a good accuracy in all orientations using uniform mesh, while an aggressive time stepping can be used. Compared with the previous CA models, the VISCA method exhibits a significant reduction in computational cost and improvement in computational efficiency. The nucleation is simulated based on a continuous nucleation model, and the Kurz-Giovanola-Trivedi (KGT) model is used to simulate the growth of the grains nucleated in the bulk material. The proposed VISCA method is comprehensively validated using an analytical solution, a casting simulation, and a directional solidification experiment. A quantitative comparison with the experiment in the directional solidification example shows that the difference in position of Columnar-to-Equiaxed Transition (CET) between the simulation and the experiment is within 1 mm while the length of directional solidification is 60 mm.

1. Introduction

The microstructure of an alloy has a significant influence on its mechanical and functional properties. The microstructure varies according to different manufacturing processes such as welding, casting, and additive manufacturing. A review of the solidification microstructure by Boettinger et al.¹⁾ describes that the material microstructure is a strategic link among its processing, properties, and performance. Thus, the microstructure simulation becomes the critical point in the prediction of mechanical and functional properties of the produced parts. In addition to experiments, the dynamic microstructural evolution can be obtained through numerical simulation with the development of computing technology.

As summarized by Zinovieva et al.,²⁾ the most popular methods for the simulation of microstructure evolution during solidification are phase-field (PF) and cellular automata (CA). The PF method³⁾ is an approach that can solve different kinds of problems about phase interface evolution. For the solidification simulation, the PF method performs well with the reproduction and prediction of the microstructure evolution at the microscale.⁴⁾ Based on the work of kobayashi,⁵⁾ Karma proposed a quantitative Phase field model for pure materials^6,7) and binary alloys^8,9) by thin interface analysis. The computational efficiency extremely increased compared to the previous models due to the extended interface thickness. However, due to the enormous demands of computational effort within a mesoscale domain (several dendrites), it is difficult to simulate the microstructural evolution at a larger scale, such as a weld pool or an ingot.

The CA method is known as another well-developed mathematical method for simulating metal solidification. It is commonly used in macroscale and mesoscale solidification simulations. In mesoscale, Beltran-Sanchez and Stefanescu^10,11) proposed a CA model without using the KGT (Kurz et al.¹²⁾ in 1986) model. They introduced the sharp interface model to control the phase boundary instead of integrating the tip velocity to compute the length of the primary dendrite. Then, Zhu and Stefanescu¹³⁾ compared the Lipton-Glicksman-Kurz (LGK) model (Lipton et al.¹⁴⁾ in 1984) with the virtual front tracking method, which also introduced the sharp interface model in the CA method. The virtual front tracking method shows a good agreement with the LGK model on the dendritic morphology in 2D and 3D simulations. Zhang and Sun¹⁵⁾ coupled the same CA method with the Lattice Boltzmann method to simulate the microporosity formation and dendrite growth during solidification. Wei¹⁶⁾ introduced a height function method to improve the calculation accurucy of interface curvature. Jelinek and Felicelli¹⁷⁾ coupled the CA method with the Lattice Boltzmann method to simulate dendritic solidification under forced convection and the computer program was parallelized using the Message Passing Interface (MPI) technique. These CA models mainly focus on mesoscale simulations, and they are challenging to be extended to larger scales.

On a larger scale, an original CA method was first proposed in 1993 to reproduce the different preferred growth orientations during nucleation in an isothermal environment by Rappaz and Gandin.¹⁸⁾ It was realized that the squared mesh would misdirect the primary dendrite orientation resulting in artificial anisotropy, and it is necessary to correct this misdirection after a few computation steps. The inherent initiator of anisotropy is induced by the neighboring cell capture method coupled with the computation mesh. The spreading speed is anisotropic in Von Neumann and Moore neighboring cells, and they will twist the primary dendrite orientation to 0° and 45°, respectively. More details of this anisotropy can be found in the article of Marek.¹⁹⁾ In addition, to address the situation of a non-isothermal environment, Rappaz et al.²⁰⁾ used a finite element (FE) method to compute the temperature field in the domain. The FE method splits the domain into elements, and there are multiple CA cells located in every single element. Because of the coarser mesh density of the heat transfer simulation, larger scales can be reproduced. Their CA–FE model deeply influenced the solidification simulation in the following decades. Ishida et al.²¹⁾ used Rappaz’s model investigate and optimum pouring temperature of the molten for obtaining the fine grain structure during high Mn steel casting. Wang et al.²²⁾ investigated columnar to equiaxed transition in free-cutting steel based on the CA–FE model. Zhu et al.²³⁾ investigated the effect of solidification pressure on the compactness degree of high nitrogen steel with the CA–FE model. Zhu et al.²⁴⁾ proposed a modified CA method in 2001, in which the temperature field and the concentration field were computed using a finite-difference (FD) method. Zhu considered the curvature undercooling by computing the local solid fraction at the solid-liquid interface. The curvature undercooling is considered in the computation of tip dendrite growth velocity using the KGT model. However, in their work, the correction of the mesh anisotropy was still not considered, and it twisted all the dendrites to 0° or 45°.

In recent years, various researchers have contributed to developing the CA method proposed by Rappaz¹⁸⁾ to reduce anisotropy. For instance, Luo et al.²⁵⁾ developed the CA method by proposed a modified decentered square growth algorithm to reduce the mesh anisotropy. Zhan et al.²⁶⁾ reported a multilayer mesh method to simulate the dendritic growth with different orientations. This approach is highly time-consuming with complicated operations due to the usage of multilayer mesh. To improve the computational efficiency, Zinovieva et al.²⁾ designed a modified technique to realize the correction in a uniform mesh of squared cells. They applied the modified method for predicting the grain structure evolution of 316L austenitic stainless steel under selective laser melting process.²⁷⁾ The simulation results are in good qualitative agreement with experimental data. However, the distance to trigger the next cell is not constant in their algorithm, and some cells in specific orientations will be triggered at a longer distance. The computational accuracy of the CA method is relative to the maximum value of trigger distance among all different orientations, while the computational efficiency is relative to its minimum value. As a result, the balance between accuracy and efficiency is a bottleneck to the simulations in reducing anisotropy, while scarce attention has been paid to this aspect.

Based on the CA model proposed by Rappaz and Gandin,^18,20) a high-efficiency virtual submesh cellular automata (VISCA) method is developed in this work. The uniform mesh of squared cells is used with the same computational accuracy for each growth direction. The CA mesh is conceptually divided into a series of submeshes by giving a unique index to each CA cell. The continuous nucleation model and the KGT model are employed for the simulation of nucleation and crystal growth.

The paper is outlined as follows. Section 2 introduces the heat transfer, nucleation and grain growth formulations used in the VISCA method. The method of generating submeshes is explained in this section as well. Section 3 describes the numerical implementation of the VISCA method and illustrates the neighboring cell capture method during grain growth in detail. Section 4 validates this model using analytical, numerical and experimental results from different solidification cases, and the substantial efficiency improvement is also illustrated in this section. Section 5 draws the conclusion and points out future research directions.

2. Model Description

The solidification process consists of nucleation and grain growth. According to the previous work by Rappaz and Gandin,¹⁸⁾ there are three determining factors for the final solidified microstructure: (i) the heterogeneous nucleation on the surface and within the molten metal, (ii) the crystallographic orientation generated randomly during nucleation, and (iii) the dendrite tip growth velocity, which is determined by local undercooling. The two-dimensional solidification process is schematically illustrated in Fig. 1, where θ is the angle between the preferred crystallographic direction and the global x-axis. The dotted line represents the solid/liquid interface.

Fig. 1.

Grains growth with different orientations. (Online version in color.)

2.1. Virtual Submesh

Because of the mesh anisotropy caused by squared cells, the simulation of grains growth at different directions is unrealizable in a single uniform mesh. In this work, as is shown in Fig. 2(a), the global mesh is conceptually divided into a series of virtual submeshes with different orientations, and the grain is located within the submesh with an orientation parallel to it. For the simulation of a crystal that always grows along with <100> orientation (such as nickel-base superalloy and aluminum-silicon alloy), the growth anisotropy (not the mesh anisotropy) is fourfold symmetrical in the two-dimensional case, and the nucleation orientation can be expressed within (0, π/2). As shown in Fig. 2(a) (i)–(vi), six classes of virtual submeshes are able to realize 24 classes of nucleation orientation within (0, 2π). After the nucleation calculation, a random number θ indicating the crystal orientation degree continuously that ranges from 0 to π/2 is given to a series of cells (nucleation points) in the liquid. A random number is sorted into six types of virtual submeshes by a classification rule shown in Fig. 2(b). The submesh is then generated simultaneously around the nucleation point, and more details are provided in Section 3.

Fig. 2.

(a) Virtual submeshes for different orientations and (b) its classification rule. The ratio between the virtual submesh step length of ΔX and the global mesh step length of Δx is shown under the submesh in Fig. 2(a). The continuous orientation degree θ is sorted into six classes of submesh by a classification rule shown in Fig. 2(b).

The virtual submesh step length ΔX is larger than the global mesh step length Δx, and the ratio between them is defined as ϵ. The ϵ of each virtual submesh is independent, and the maximum and minimum values are defined as ϵ_max and ϵ_min, respectively. The difference between ϵ_max and ϵ_min is ϵ_delta, and the average value of ϵ is ϵ_mean. On one hand, for a constant Δx, a larger ϵ_max will obviously decrease the computational accuracy due to fewer cells in the same space. On the other hand, a smaller ϵ_min will intensively increase the lower limit of the time step, which can decrease the computational efficiency. When the ϵ_delta equals zero (the computational accuracy is the same at different orientations), the resultant computational efficiency will be higher. In the work of Zinovieva²⁾ and Zhan,²⁶⁾ the trigger distance between cells in a direction can be several times greater than that in the other directions. In other words, the computational accuracy in their algorithm is uneven for every grain orientation. Consequently, the lower limit of the time step has to be set to a lower level, which is already unnecessary in some directions.

The CA method can simulate with a limited number of orientations n_2π. A larger value of n_2π will increase the magnitude of ϵ_mean. Furthermore, a smaller ϵ_delta will also increase the magnitude of ϵ_mean. Obviously, the size of the global matrix is positively correlated with ϵ_mean, and a more prominent global matrix will cost more CPU and memory resources. To maintain a balance of reality and computational efficiency, the number of orientations n_2π should be more significant, and the values of ϵ_delta and ϵ_mean should be smaller. Because the 24 classes of orientation are sufficient to reproduce the real solidification process, the n_2π is set to 24. The value of ϵ_mean is set to approximately 6 to ensure that the ϵ_delta is always smaller than one under the current n_2π in the VISCA method. The precise value of ϵ at every orientation is shown in Fig. 2(a).

2.2. Heat Transfer Formulation

In this work, the temperature field is solved using a finite difference method. The governing equation for the heat transfer process can be expressed as follows:   

$$\rho C_p\frac{\partial T}{\partial t}=\nabla \cdot (\lambda \nabla T)+\rho L\frac{\partial f_s}{\partial t}+\varphi ,$$(1)

where T is the temperature, t is the time, λ is the thermal conductivity, ρ is the density, C_p is the heat capacity, L is latent heat, f_s is solid fraction, ϕ is the heat source term. Material properties are interpolated using f_s between solid and liquid phases. The boundary condition is specified as T_b = T^evn, where T^evn is the environment temperature, and T_b is the boundary temperature. Both explicit and implicit time stepping schemes¹ can be employed for solving Eq. (1). A coarser mesh can be used to solve heat equation for reducing computation cost, and intergrid transfer for interface computed from CA (by averaging) and temperature field computed from heat equation (by interpolation) is required. In this work, the grid size for solving heat equation is set to be ten times as large as that of the virtual submesh.

¹  To improve the computation stability of explicit schemes, the finite difference form in the enthalpy method by Voller²⁸⁾ can be used for computing derivatives in Eq. (1).

2.3. Nucleation Formulation

Following the pragmatic approach used by Rappaz et al.,¹⁸⁾ the continuous nucleation distribution can be expressed as follows:   

$$n(\mathrm{\Delta }T)={\int }_0^{\mathrm{\Delta }T}\frac{n_{max}}{\mathrm{\Delta }{T}_{\sigma }\sqrt{2\pi }}\cdot e^{\left(-\frac{{(\mathrm{\Delta }T-\mathrm{\Delta }{T}_N\hspace{0.17em})\hspace{1em}}^2}{2\mathrm{\Delta }{T}_{\sigma }{}^2}\right)}\mathrm{d}(\mathrm{\Delta }T),$$(2)

where n is the nucleation density, ΔT is the undercooling, n_max is the maximum nucleation density, ΔT_N is the mean undercooling, and ΔT_σ is the standard deviation, and the maximum nucleation density in 2D can be derived from the work by Zinovieva²⁾ as n_max = $\sqrt{\frac{\pi }{6}}\cdot {(n_{max}{}_{2D})}^{\frac{3}{2}}$. The nucleation density increment between tⁿ and tⁿ⁺¹ time step can be expressed as Δn = n(ΔTⁿ⁺¹) − n(ΔTⁿ). For a nonisothermal situation, the nucleation distribution is not uniform in space. By solving the Eq. (2), the nucleation probability function (2D version) is constructed for every global mesh cell, which is given as follows:   

$$\begin{array}{l}{P}_n=\mathrm{\Delta }n\cdot S_{CA}=\\ \frac{n_{max2D}\hspace{1em}\cdot \mathrm{\Delta }{x}^2}{2}\left(erf\left(\frac{\mathrm{\Delta }{T}_N-\mathrm{\Delta }{T}^n}{\sqrt{2}\mathrm{\Delta }{T}_{\sigma }}\right)-erf\left(\frac{\mathrm{\Delta }{T}_N-\mathrm{\Delta }{T}^{n+1}}{\sqrt{2}\mathrm{\Delta }{T}_{\sigma }}\right)\right),\end{array}$$(3)

where S_CA is the area of a global mesh equal to Δx², and P_n is the nucleation probability for correlative cells at the current time step. To reduce the number of random numbers, only correlative cells (0 ≤ ΔTⁿ < ΔTⁿ⁺¹ ≤ (ΔT_N + ΔT_σ)) are activated for nucleation.

2.4. Grain Growth

The growth velocity of the primary dendrite tip V_tip can be considered as a function of local undercooling ΔT. The equations from the LGK model¹⁴⁾ and KGT model¹²⁾ are used to gain the relationship between ΔT and V_tip. The following equations¹²⁾ link the V_tip and the Peclet number P_c:   

$$V_{tip}{}^2\cdot \frac{{\pi }^2\cdot \mathrm{\Gamma }}{P_c^2\cdot D_l^2}+V_{tip}\cdot \frac{m_l\cdot C_0\cdot (1-k)\cdot {\xi }_c}{D{}_l\cdot [1-(1-k)\cdot Iv\text{(}{P}_c\text{)}]}+G=0$$(4)

  

$${\xi }_c=1-\frac{2k}{\sqrt{\text{1}+{\left(\frac{\text{2}\pi }{P_c}\right)}^2}-1+2k},$$(5)

where Γ is the Gibbs–Thomson coefficient, D_l is the liquid diffusion coefficient, k is the partition coefficient, m_l is the liquidus slope, C₀ is the initial concentration, G is the temperature gradient, P_c is the concentration Peclet number, and Iv(P_c) is the Ivanstov function of concentration Peclet number expressed as follows:   

$$Iv\text{(}{P}_c\text{)}=P_c\cdot e^{P_c}\cdot {\int }_{P_c}^{\infty }\frac{e^{-x}}{x}\hspace{0.17em}\text{d}x\text{.}$$(6)

By solving Eq. (2) with a given value of V_tip, the concentration Peclet number P_c can be calculated. The following equations¹⁴⁾ link the P_c and the undercooling ΔT:   

$$\mathrm{\Delta }T=m_l\cdot C_{\text{0}}\cdot \left(\text{1}-\frac{1}{\text{1}-\text{(1}-k\text{)}\cdot Iv\text{(}{P}_c\text{) }}\right)+\frac{\text{2}\mathrm{\Gamma }}{R}$$(7)

  

$$P_c=\frac{V_{tip}\cdot R}{\text{2}{D}_l},$$(8)

where R is the dendrite tip radius and Eq. (7) is satisfied only when the liquidus slope and partition coefficient are constant. The undercooling ΔT is solved by substituting P_c (which is calculated above) into Eqs. (7) and (8). By solving these equations with parameters in Table 2, the relationship between V_tip and ΔT is fitted in Eq. (9) and plotted in Fig. 3.   

$$V_{tip}=\text{2}\text{.499}\times {\text{10}}^{-\text{5}}\mathrm{\Delta }T+\text{1}\text{.551}\times {\text{10}}^{-\text{6}}\mathrm{\Delta }{T}^{\text{2}}+\text{5}\text{.606}\times {\text{10}}^{-\text{7}}\mathrm{\Delta }{T}^{\text{3}}$$(9)

Table 2. Material properties of Al−7mass%Si used in the simulation.

Property	Variable	Unit	Value
Density	ρ	kg/m3	2600
Heat capacity	Cp	J/(kg · K)	998
Solid thermal conductivity	λs	W/(m · K)	170
Liquid thermal conductivity	λl	W/(m · K)	66
Liquidus temperature	Tl	K	930
Mean undercooling	ΔTN	K	8
Standard deviation	ΔTσ	K	4
Maximum nucleation density	nmax	1/m3	9 × 1010
Liquidus slope	ml	K/wt%	−6
Initial concentration	C0	wt%	7
Partition coefficient	k	1	0.117
Gibbs-Thomson coefficient	Γ	K · m	9 × 10−8
Diffusion coefficient	Dl	m2/s	1 × 10−9

Fig. 3.

The relationship between V_tip and ΔT is fitted by 1000 couples of discrete solutions. The value of V_tip ranges from 1 um/s to 0.1 m/s and densified around 1 um/s. (Online version in color.)

2.5. Time Increment

There are two kinds of time increment used during the entire calculation. The CA time increment Δt is as follows:   

$$\mathrm{\Delta }t=\text{0}\text{.25}\cdot \text{min}\left(\frac{\mathrm{\Delta }X}{\text{max(}{V}_{tip}\text{)}},\frac{\mathrm{\Delta }{X}^2}{D_l}\right),$$(10)

where max(V_tip) is the maximum value of the dendrite tip velocity of the entire computational domain. The global time increment Δt is the duration of one complete calculation cycle, including one cycle of temperature, nucleation, and growth simulations. The time increment for heat transfer is denoted as Δt_T, and for explicit time stepping Δt_T can be expressed as follows:   

$$\mathrm{\Delta }{t}_T=\frac{\mathrm{\Delta }{x}_T{}^{\text{2}}}{\text{4}\frac{\lambda }{\rho C_p}}=\frac{\mathrm{\Delta }t}{r},$$(11)

where Δx_T is the mesh step length for heat equation. Since a smaller time step has to be employed for explicit heat equation, the temperature field is solved r times during Δt to ensure synchronization in time. For the implicit time stepping, Eq. (10) is used, and the Δt_T is equal to Δt.

3. Procedure of the Virtual Submesh Method

The computation procedure is illustrated in Fig. 4.

Fig. 4.

Calculation procedure of solidification simulation of the VISCA method. The calculation procedure consists of three parts: heat diffusion calculation, nucleation calculation, and grain growth calculation.

The first part is the simulation of the temperature field using an updated interface (i.e., values of f_s) computed from the previous time step. The second part of the simulation procedure is the nucleation calculation. The nucleation probability of each cell in the global mesh can be derived from Eq. (3) through the local undercooling ΔT. After rolling a random number, the new nucleation cell in the liquid at the current time step is selected with a random orientation θ. Then, a virtual submesh with the same orientation is constructed around the nucleated cell. The third part of the simulation procedure is the grain growth calculation. The virtual submesh around the nucleation cell expands along with the grain growth, as shown in Fig. 5, where L_d is the primary dendrite length of the current cell (not the real one), and ΔX is the virtual submesh step length. After the nucleation, the growth state L_d is stored in a specific matrix. At the beginning of the third part, the nucleation points on the global mesh will be led to three different ways according to the specific conditions:

Fig. 5.

Schematic diagrams of VISCA growth and the generation of virtual submesh: (a) When L_d of the central cell is smaller than ΔX, a₁ a₄ are still not activated (b) When L_d of the central cell is more significant than ΔX, a₁a₄ are activated, and the L_d of them starts to grow. The white cells are global meshes, and the black ones are virtual submeshes. The gray ones are skin cells of the growing grain.

(a) If the L_d/ΔX of the current cell is within the range of (0, 1), the V_tip can be obtained from Eq. (9). The growth of dendrite length at the current time step will be added to L_d (as shown in Fig. 5(a)), which is the product between V_tip and Δt.

(b) For L_d/ΔX larger than 1, the L_d of the current cell will be set to ΔX, and the adjacent cells of the inactivated state in virtual submesh (a₁–a₄ in Fig. 5(a)) will be activated, and the surplus of L_d − ΔX will be stored in the growth state L_d matrix. Afterwards, a series of cells on the global matrix (the offset address of them are shown in Table 1) surrounded by four adjacent ones on the submesh will be given the same value as the central cell (the value is orientation degree θ), which can be named as the cover skin (because the submesh is not discontinuous on the global mesh, the grain should be performed on global mesh). Then, the virtual submesh expands around the activated cells. The offset address for submesh expansion is shown in Table 1, Fig. 5(b).

Table 1. Offset address for virtual submesh expansion.

Submesh index	Orientation degree	Offset address
i	(0, π/24) ∪ (11π/24, π/2)	(6,0) (0,6) (−6,0) (0,−6)
ii	(π/24, π/8)	(6,2) (−2,6) (−6,−2) (2,−6)
iii	(π/8, 5π/24)	(5,3) (−3,5) (−5,−3) (3,−5)
iv	(5π/24, 7π/24)	(4,4) (−4,4) (−4,−4) (4,−4)
v	(7π/24, 3π/8)	(3,5) (−5,3) (−3,−5) (5,−3)
vi	(3π/8, 11π/24)	(2,6) (−6,2) (−2,−6) (6,−2)

(c) Once the submesh expanding into skin cells, the submesh stops expanding, and the cell will not trigger any adjacent cells. As shown in Fig. 6(a), the red grains are surrounded by the submesh cells a₁–a₈, and a₁ is found to expand into the skin cells of another grain. The submesh cells b₁–b₇ in Fig. 6(b) are generated from a₁–a₈, and the submesh cell a₁ will not trigger other cells. As shown in Fig. 6(c), the cell a₁ still can cover the red skin cells except for the blue cells.

Fig. 6.

Schematic diagrams of grain growth simulation between grains: (a)–(c) showing the submesh expansion process when the growing grain (in red) touches other grains (in blue or green). (Online version in color.)

After three parts of the calculation, the solid fraction f_s will be updated, which then determines the material properties and latent heat for the next time step.

4. Results and Discussion

4.1. Validation Examples

A numerical simulation of microstructure evolution during solidification (shown in Fig. 7) is produced by the VISCA method. The values next to the colorbar are the index of phase and orientation (−1 is for liquid and (0–π/2) is for solid). Because of the abundant successful simulation experiences and the <100> growth direction, the Al−7mass%Si alloy is used in this work for the simulation of grain growth.^18,29) The material parameters used in the current simulation are shown in Table 2. The parameters of nucleation, diffusion and interfacial energy refer to the work by Zhu et al.²⁹⁾ The thermophysical parameters refer to the work by Liu et al.^30,31)

Fig. 7.

Evolution of microstructure computed by the VISCA algorithm at different stages. The values of the colorbar are the index of phase (−1 is for liquid and (0–π/2) is for solid of different orientations). (Online version in color.)

The metal melt of Al−7mass%Si freeze in an atmosphere of 300 K, and the initial temperature of the melt is 1000 K. The computational domain is two-dimensional with a size of 10 mm × 10 mm. The global CA mesh size is 1000 × 1000 and the temperature mesh size is 100 × 100. The implicit algorithm is chosen for this case. With the rapid decrease in melt temperature, the nucleation and grain growth process originate almost simultaneously. The topography of crystal grains is visually similar to the cast experiment.³²⁾

To further validate the VISCA algorithm, a simulation of directional solidification is implemented for comparison with the experiment and simulation published by Liu et al.³⁰⁾ The experiment is implemented in a low gravity environment where natural convection does not affect the nucleation and temperature field, and it is more suitable for validating the solidification simulation algorithm. The ‘FM1’ experiment in their work is simulated to verify the VISCA algorithm. The directional solidification experiment includes three cooling stages with different pulling velocity: 0.01 mm/s in stage I, 0.2 mm/s in stage II and 3 mm/s in stage III. The thermophysical and nucleation parameters used in the simulation are the same as the simulation work of Liu. The computational domain is two-dimensional with a size of 4 mm × 60 mm, and the global mesh size is 200 × 3000. The vertical temperature field used in simulation is generated from the cooling curves measured at several evenly distributed locations (positions of TC6–TC12 are plotted in Fig. 8(a)) during experiment, and the horizontal temperature field is considered to be uniform. Some points on the cooling curves (endpoints of stage I–stage III) are taken to generate simplified cooling curves by linear interpolation (plotted in Fig. 8(b)).

Fig. 8.

The FM1 experiment (in the work of Liu) is simulated to verify the VISCA algorithm: (a) The directional solidification is simulated by the VISCA algorithm. The positions of thermocouples (TC6–TC12) are plotted and the position of Columnar-to-Equiaxed Transition is at 127 mm, which is the same as the experiment. (b) The cooling curves are generated by thermocouples. (Online version in color.)

The simulation result is shown in Fig. 8(a) with a length scale and a colorbar for grain orientation (range from 0 to π/2). The position of Columnar-to-Equiaxed Transition (CET) is at 127 mm, which is the same as the experiment. The position of CET in simulation fits well with the experiment, and the grain growth velocity spans from 3.5 um/s to 200 um/s in this case, proving that the algorithm is able to reproduce the solidification process at both high and low growth velocities.

4.2. Reduction of Artificial Anisotropy

To assess the performance of reducing artificial anisotropy using the VISCA algorithm, the grains growth in an isothermal environment is simulated. The heterogeneous nuclei of six orientation angles θ₁–θ₆ (0, π/12, π/6, π/4, π/3, 5π/12) are uniformly located in the domain. The same case also implemented by ‘correction 2’ model of Zinovieva et al.²⁾ The dimensionless computational parameters are shown in Table 3. To avoid the loss of the primary dendrite length, The time increment Δt of VISCA is set to maximum allowed value of $\frac{\text{0}\text{.25}\mathrm{\Delta }X}{\text{max(}{V}_{tip}\text{)}}$ = 1.5, and Δt of previous method is set to maximum allowed value of $\frac{\text{0}\text{.25}\mathrm{\Delta }x}{\text{max(}{V}_{tip}\text{)}}$ = 0.25. The primary dendrite lengths can be measured by calculating the area of grains before any possible collision. The analytical solution of primary dendrite length is y = x since all the grains are grown with a constant velocity 1. For VISCA algorithm with Δt = 1.5 in Figs. 9(a), 9(b), it is obvious that the dendrite lengths of different grains remain the same and all the curves converge to the analytical solution. The results show that the velocity loss of growth is avoided. For the previous method with Δt = 0.25 in Fig. 9(c), the simulation results are similar to VISCA. The primary dendrite lengths in Fig. 9(d) also converge to the analytical solution. But the vibration frequencies of these curves are different, and it’s the reason why previous method must work with a smaller time increment. These numerical experiments illustrate that both two methods overcome artificial anisotropy, and the VISCA algorithm is able to work in a much bigger time increment Δt.

Table 3. Offset address for virtual submesh expansion.

Computational parameter	Variable	Value	Unit
Global mesh step length	Δx	1	1
Virtual submesh step length	ΔX	f (θ) ≈ 6Δx	1
Global mesh size	m × n	900 × 1200	1
Dendrite tip velocity	Vtip	1	1
Time increment (VISCA)	Δt	$0.25\frac{\mathrm{\Delta }X}{\text{max(}{V}_{tip}\text{)}}=1.5$	1
Time increment (Previous)	Δt	$0.25\frac{\mathrm{\Delta }x}{\text{max(}{V}_{tip}\text{)}}=\text{0}\text{.25}$	1

Fig. 9.

Comparison of grain growth velocity difference among grains of different orientations between VISCA and previous algorithm: (a)–(b) The primary dendrite length of growing grains by VISCA algorithm with Δt = 1.5. (c)–(d) The primary dendrite length of growing grains by previous algorithm with Δt = 0.25. The heterogeneous nuclei of six orientation angles θ₁−θ₆(0, π/12, π/6, π/4, π/3, 5π/12) are uniformly located in the domain. (Online version in color.)

4.3. Computational Efficiency

In this section, the comparison of computational efficiency between the VISCA algorithm and the ‘correction 2’ model of Zinovieva²⁾ is performed. The case used for comparison is the same as Section 4.2. The computational efficiency comparison is plotted in Fig. 10. The relationship between physical time and computation time is shown in Fig. 10(a). It is evident that the VISCA method can calculate longer physical time during the same computation time. In other words, the VISCA method has higher computational efficiency. The relative computational efficiency plotted in Fig. 10(b) is the ratio between the physical time of the VISCA algorithm and that of the previous algorithm under the same CPU time. The computational efficiency is dramatically increased by over 3.5 times using the VISCA algorithm in this work, and the computational efficiency can be simultaneously improved along with the grain growth. The efficiency advantage of the VISCA algorithm is more remarkable for the initially nucleated grains, which is approximately four times higher than that reported in previous work.²⁾ The efficiency advantage will be continually improved with grain growth.

Fig. 10.

Computational efficiency comparison between the VISCA algorithm and the previous algorithm: (a) The relationship between physical time and CPU time of each algorithm. (b) The ratio between the computational efficiency of the current and previous algorithm. (Online version in color.)

As stated above, the simulation requires correction to eliminate the velocity loss. The practical methods are increasing the mesh step length or reducing the time step. Although the chosen time step has a significant effect on the computation time, reducing the time step is still the only option to maintain the computational accuracy in the current scale. The final time step Δt is 1.5 (dimensionless) and 0.25 (dimensionless), respectively, to eliminate the grain size inhomogeneity at different orientations. In the previous algorithm, the computation accuracy at each direction is different since the distinct trigger distance between cells. For example, the trigger distance for θ equals to arctan(3/5) is $\sqrt{34}$ times higher than that for θ of 0 in the previous algorithm. A smaller trigger distance can lead to more iterations. As a result, the direction with lower accuracy always waits for the direction with higher accuracy. It will waste lots of computational effort because the total accuracy equals the lowest.

The VISCA algorithm divides the global mesh into a series of submesh with approximately equal submesh step lengths, and thus, the upper limits of the theoretical time step for grains of different orientations are similar in value. Compared to the previous algorithm using irregular accuracy for different growth directions, the submesh method will save much computational effort for those grains of excessive computational accuracy when other grains are growing in low accuracy.

5. Conclusions

Based on the cellular automata method proposed by Rappaz,¹⁸⁾ the current work provides a novel algorithm that can compute grain growth at different orientations without any direction correction. Using a virtually divided global mesh, a more aggressive time stepping can be achieved while the computational accuracy for different directions is retained the same. In traditional approaches, the computational efforts substantially increase as grains grow, eventually resulting in extremely slow computation. However, the proposed VISCA method solves this issue by reducing massive redundant time steps in specific directions and tremendously improves the overall efficiency. The simulation results from the VISCA method present good agreements with other analytical, simulation and experiment results. We anticipate deploying this framework to solidification simulations in steel casting and additive manufacturing in the future.

Acknowledgments

This work was jointly funded by the National Science and Technology Major Project “Aeroengine and Gas Turbine” of China (2017-VII-0008-0102), the “Sailing Program” of Shanghai Science and Technology Committee of China (19YF1415900), the National Natural Science Foundation of China (Nos. 52001191), the Shanghai Municipal Science and Technology Commission Grant of China (No. 20511107700, No. 21511103600), and Shanghai Rising-Star Program (20QA1403800).

References

• 1)   W. J.  Boettinger,  S. R.  Coriell,  A. L.  Greer,  A.  Karma,  W.  Kurz,  M.  Rappaz and  R.  Trivedi: Acta Mater., 48 (2000), 43. https://doi.org/10.1016/S1359-6454(99)00287-6
• 2)   O.  Zinovieva,  A.  Zinoviev,  V.  Ploshikhin,  V.  Romanova and  R.  Balokhonov: Comput. Mater. Sci., 108 (2015), 168. https://doi.org/10.1016/j.commatsci.2015.06.026
• 3)   I.  Steinbach,  F.  Pezzolla,  B.  Nestler,  M.  Seeßelberg,  R.  Prieler,  G. J.  Schmitz and  J. L. L.  Rezende: Physica D, 94 (1996), 135. https://doi.org/10.1016/0167-2789(95)00298-7
• 4)   W. J.  Boettinger,  J. A.  Warren,  C.  Beckermann and  A.  Karma: Annu. Rev. Mater. Res., 32 (2002), 163. https://doi.org/10.1146/annurev.matsci.32.101901.155803
• 5)   R.  Kobayashi: Physica D, 63 (1993), 410. https://doi.org/10.1016/0167-2789(93)90120-P
• 6)   A.  Karma and  W.-J.  Rappel: Phys. Rev. E, 53 (1996), R3017. https://doi.org/10.1103/PhysRevE.53.R3017
• 7)   A.  Karma and  W.-J.  Rappel: Phys. Rev. E, 57 (1998), 4323. https://doi.org/10.1103/PhysRevE.57.4323
• 8)   A.  Karma: Phys. Rev. Lett., 87 (2001), 115701. https://doi.org/10.1103/PhysRevLett.87.115701
• 9)   B.  Echebarria,  R.  Folch,  A.  Karma and  M.  Plapp: Phys. Rev. E, 70 (2004), 061604. https://doi.org/10.1103/PhysRevE.70.061604
• 10)   L.  Beltran-Sanchez and  D. M.  Stefanescu: Metall. Mater. Trans. A, 34 (2003), 367. https://doi.org/10.1007/s11661-003-0338-z
• 11)   L.  Beltran-Sanchez and  D. M.  Stefanescu: Metall. Mater. Trans. A, 35 (2004), 2471. https://doi.org/10.1007/s11661-006-0227-3
• 12)   W.  Kurz,  B.  Giovanola and  R.  Trivedi: Acta Metall., 34 (1986), 823. https://doi.org/10.1016/0001-6160(86)90056-8
• 13)   M.  Zhu and  D. M.  Stefanescu: Acta Mater., 55 (2007), 1741. https://doi.org/10.1016/j.actamat.2006.10.037
• 14)   J.  Lipton,  M. E.  Glicksman and  W.  Kurz: Mater. Sci. Eng., 65 (1984), 57. https://doi.org/10.1016/0025-5416(84)90199-X
• 15)   Q.  Zhang,  D.  Sun,  S.  Pan and  M.  Zhu: Int. J. Heat Mass Transf., 146 (2020), 118838. https://doi.org/10.1016/j.ijheatmasstransfer.2019.118838
• 16)   L.  Wei,  Y.  Cao,  X.  Lin,  M.  Wang and  W.  Huang: Comput. Mater. Sci., 156 (2019), 157. https://doi.org/10.1016/j.commatsci.2018.09.019
• 17)   B.  Jelinek,  M.  Eshraghi,  S.  Felicelli and  J. F.  Peters: Comput. Phys. Commun., 185 (2014), 939. https://doi.org/10.1016/j.cpc.2013.09.013
• 18)   M.  Rappaz and  C. A.  Gandin: Acta Metall. Mater., 41 (1993), 345. https://doi.org/10.1016/0956-7151(93)90065-Z
• 19)   M.  Marek: Physica D, 253 (2013), 73. https://doi.org/10.1016/j.physd.2013.03.005
• 20)   C. A.  Gandin and  M.  Rappaz: Acta Metall. Mater., 42 (1994), 2233. https://doi.org/10.1016/0956-7151(94)90302-6
• 21)   H.  Ishida,  Y.  Natsume and  K.  Ohsasa: ISIJ Int., 48 (2008), 1728. https://doi.org/10.2355/isijinternational.48.1728
• 22)   J.  Wang,  F.  Wang,  C.  Li and  J.  Zhang: Steel Res. Int., 81 (2010), 150. https://doi.org/10.1002/srin.200900109
• 23)   H.-C.  Zhu,  Z.-H.  Jiang,  H.-B.  Li,  J.-H.  Zhu,  H.  Feng,  S.-C.  Zhang,  B.-B.  Zhang,  P.-B.  Wang and  G.-H.  Liu: Steel Res. Int., 88 (2017), 1600509. https://doi.org/10.1002/srin.201600509
• 24)   M.  Zhu and  C.  Hong: ISIJ Int., 41 (2001), 436. https://doi.org/10.2355/isijinternational.41.436
• 25)   S.  Luo and  M.  Zhu: Comput. Mater. Sci., 71 (2013), 10. https://doi.org/10.1016/j.commatsci.2012.12.040
• 26)   X.  Zhan,  Y.  Wei and  Z.  Dong: J. Mater. Process. Technol., 208 (2008), 1. https://doi.org/10.1016/j.jmatprotec.2007.12.130
• 27)   A.  Zinoviev,  O.  Zinovieva,  V.  Ploshikhin,  V.  Romanova and  R.  Balokhonov: Mater. Des., 106 (2016), 321. https://doi.org/10.1016/j.matdes.2016.05.125
• 28)   V. R.  Voller: Int. J. Heat Mass Transf., 51 (2008), 823. https://doi.org/10.1016/j.ijheatmasstransfer.2007.04.025
• 29)   M.  Zhu,  J. M.  Kim and  C.  Hong: ISIJ Int., 41 (2001), 992. https://doi.org/10.2355/isijinternational.41.992
• 30)   D.  Liu,  N.  Mangelinck-Noël,  C. A.  Gandin,  G.  Zimmermann,  L.  Sturz,  H.  Nguyen-Thi and  B.  Billia: Acta Mater., 93 (2015), 24. https://doi.org/10.1016/j.actamat.2015.03.058
• 31)   D.  Liu,  E.  Guo,  L.  Wang,  F.  Kang and  Q.  Lan: Int. J. Cast Met. Res., 20 (2007), 254. https://doi.org/10.1179/136404607X268085
• 32)   Y.  Birol: Mater. Sci. Technol., 28 (2012), 363. https://doi.org/10.1179/1743284711Y.0000000041