Numerical Simulation of Three-Dimensional Mesoscopic Grain Evolution: Model Development, Validation, and Application to Nickel-Based Superalloys

The mesoscopic grain model is a multiscale model which takes into account both the dendrite growth mechanism and the vast numerical computation of the actual castings. Due to the pursuit of efficient computation, the mesoscopic grain calculation accuracy is lower than that of dendrite growth model. Improving the accuracy of mesoscopic grain model is a problem to be solved urgently. In this study, referring to the calculation method of solid fraction in microscopic dendrite growth model, a cellular automata model of 3D mesoscopic grain evolution for solid fraction calculated quantitatively at the scale of cell is developed. The developed model and algorithm validation for grain growth simulation is made by comparing the numerical results with the benchmark experimental data and the analytical predictions. The results show that the 3D grain envelopes simulated by the developed model and algorithm are coincident with the shape predicted by the analytical model to a certain extent. Then, the developed model is applied to the numerical simulation of solidification process of nickel-based superalloys, including equiaxed and columnar dendritic grain growth. Our results show good agreement with the related literature.


Introduction
Grain evolutions including nucleation and growth during solidification process have a profound effect on the material mechanical properties of alloys. In particular, the competitive growth of equiaxed and columnar dendritic grains determine the microstructure of alloys. Over the past few decades, many computer models and methods have been developed to study theoretically and experimentally the grain growth in alloy solidification, such as the volume-averaged, front tracking, cellular automata (CA), and phase field models [1][2][3][4][5]. Among them, due to the advantages of clearly describing the physical phenomena and high computational efficiency of solidification transition, the CA model has been widely used and developed for prediction of grain growth in additive manufacturing [6,7], welding [8,9], and casting [10,11].
In 1993, the CA model was applied for the first time in alloy solidification by Rappaz and Gandin [12]. It mainly dealt with nucleation and growth of grains in a uniform temperature field. At this point, the basic framework of the CA model for grain growth simulation was constructed, although it still had some drawbacks. Subsequently, in order to satisfy the adaptability of non-uniform temperature field at the same time, the "2D rectangular algorithm" [13] and "decentered square algorithm" [14,15] coupled with finite element (FE) calculations were proposed. Similarly, the CA model could also be coupled with the finite difference (FD) method to form the CA-FD model [16], or coupled with the finite volume (FV) method to form the CA-FV model [17]. As one of the effective methods to simulate the growth of grain, the CA-FE model integrated into the casting commercial software of ProCAST have been widely used [18,19].
However, the growing octahedron in these models is no longer considered if the corresponding cell is fully surrounded by the mushy cells [15]. In this way, this CA-FE model is still a weak coupling method and overestimates the heat release of cell transformation in the mushy zone [20]. In the CA model, the solid fraction of growing cell or interface cell in the mushy zone cannot be calculated quantitatively at the scale of CA grid. And it is calculated approximately at the scale of FE grid by the Scheil's equation and the Scheil microsegregation model. Recently, Guillemot et al. [20] have employed a front tracking method proposed by Gandin [21] to deal with the coupling problem. Similar works have also been conducted by Carozzani et al. [22] and Liu et al. [23,24].
Additionally, the CA model of grain growth, ignoring the dendritic morphology of the grain itself, is generally referred to as mesoscopic model. Corresponding CA models of dendritic growth called microscopic model have been proposed and developed, typically works like Nastac's [25], Zhu's [26] and Pan's [27]. In this regard, it is a full coupling model that accounts for the time evolution of the fraction of solid predicted at the scale of CA grid in the CA model of dendritic growth. It is further observed in the CA models of two scales that both are faced with the same problems of grid anisotropy [13,28], transformation rules [14,29], and coupling methods [30]. From the point of view of the simulation method, the method for dealing with the same problem in the CA models of two scales are interrelated [16,31]. Nevertheless, there is a problem that the CA model of dendritic growth requires sufficiently higher mesh resolution than the CA model of grain growth [32].
The purpose of this paper is to achieve the full coupling simulation of grain growth so as to improve the simulation accuracy. A 3D modified cellular automaton coupling model of grain growth for solid fraction calculated quantitatively at the scale of CA gird is proposed. The relevant models and algorithms are presented and validated. Then, the developed model is applied to simulate the formation of equiaxed and columnar dendritic grains during solidification process of nickel-based superalloys.

Model Description and Numerical Algorithm
A 3D computation domain is firstly divided into a series of regular hexahedral cells. The cell spacing is equal to the corresponding mesh size. Each cell is characterized with such variables as solid fraction, grain growth rate, grain orientation, etc. All the cell states are defined as three types: liquid (f s = 0), solid (f s = 1), or interface (0 < f s <1). And the cell neighbors are characterized by the first-order neighbors (6 neighborhoods), second-order neighbors (12 neighborhoods), and third-order neighbors (8 neighborhoods). At the beginning of the numerical simulation, the cells in computation domain are initialized with given initial and boundary conditions. When the undercooling reaches or exceeds the critical nucleation undercooling, the grain begins to nucleate and grow. As the dendritic grain grows, the solid fraction of growing cell or interface cell is calculated quantitatively at the scale of CA grid. The governing equations for calculating the nucleation distribution and growth rate of dendritic grains are presented. The numerical algorithm about capture rules for new interface cells and the calculation of solid fraction at the scale of CA grid are described in detail below.

Grain Nucleation
Only heterogeneous nucleation is considered in this study, since high undercooling is not achieved due to the presence of foreign phases. For the sake of practicability, the continuous nucleation model proposed by Rappaz and Gandin et al. [12] is used to describe the grain nucleation. The total grain density n(∆T) is calculated by: where ∆T is the undercooling. The continuous nucleation distribution dn/d∆T is given by: where n max is the total density of grains, ∆T σ is the standard deviation of undercooling, and ∆T N is the mean undercooling in nucleation, respectively. Additionally, the total density of grains on the two-dimensional section (n max ) 2D can be obtained by the stereological relationships [12], (n max ) 2D = (6·n max 2 /π) 1/3 .
Unlike the numerical realization of nucleation process determined by random number method, in this present study, the nucleation position and critical undercooling are randomly assigned to the nucleation region in advance, and then the actual undercooling and the critical undercooling of each cell are compared in each time step to determine the nucleation.

Grain Growth
Grain growth rate is mainly determined by the growth behavior of dendrite tips. According to the Lipton-Glicksman-Kurz (LGK) growth model [33,34], the dendrite tips growth of alloys are affected by thermal undercooling (∆T t ), constitutional undercooling (∆T c ) and curvature undercooling (∆T R ). With the growth of dendrite tips, latent heat and solute are released from tips. At the interface of dendrite tips, heat and solute transfer lead to changes in thermal and constitutional undercooling. In addition, the growth of dendrite tips also lead to the change of curvature undercooling. At a given total undercooling, the grains will grow under the interaction of various factors until they grow up completely or encounter obstacles to stop growth.
The thermal undercooling, constitutional undercooling and curvature undercooling are defined as: where Iv(P) is the Ivantsov function of the Peclet number; P t = RV/(2α) and P c = RV/(2D) are the heat Peclet number and solutal Peclet number, respectively; ∆H is the latent heat; c p is the specific heat; m L is the liquidus slope; C 0 is the initial components concentration of alloys; k 0 is the solute partition coefficient; Γ is Gibbs-Thompson coefficient; R is the tip radius; V is the interface growth rate; α is thermal diffusion coefficient; and D is solute diffusion coefficient. The 3D Ivantsov function is expressed as: where E 1 (P) is the exponential integral function. The relationship between the tip radius and the interface growth rate can be established by combining the above equations, but it cannot be solved, and other more definite conditions are needed.
Based on the stability criterion of a dendrite tip, the tip radius is equal to the shortest wavelength λ i , which can be evaluated by: where G c is the concentration gradient; G is the temperature gradient; σ* is a stability parameter, it should be noted here that the stability parameter according to the microsolvability theory varies with the degree of anisotropy of the interface energy, and its value can be obtained by the linearized solvability theory [35]. For ease of numerical calculations, the above equations are realized and solved using Matlab. The growth rate is simplified as a polynomial function of local undercooling.

Capture Rules for New Interface Cells and Calculation of Solid Fraction at the Scale of the CA Grid
A single regular octahedron is defined as a constituent element of a grain. Each cell corresponds to a regular octahedron. The sizes of regular octahedrons are characterized by the distance (S L ) between any two adjacent vertices. Its three orthogonal axes ([100]/[010]/[001]) represent the directions of dendrite tips. Meanwhile, in order to avoid the infinite growth problem of the single regular octahedron, the distance (S L ) is limited to 2·∆x (∆x is the cell size). The grain envelope, which consists of a number of regular octahedrons, is also shown as a larger regular octahedron in the uniform temperature field. A single grain whose crystallographic orientation is assumed to be the Euler angles (ϕ, θ, ψ) with respect to the axis of CA, as shown in Figure 1.
where Gc is the concentration gradient; G is the temperature gradient; σ* is a stability parameter, it should be noted here that the stability parameter according to the microsolvability theory varies with the degree of anisotropy of the interface energy, and its value can be obtained by the linearized solvability theory [35].
For ease of numerical calculations, the above equations are realized and solved using Matlab. The growth rate is simplified as a polynomial function of local undercooling.

Capture Rules for New Interface Cells and Calculation of Solid Fraction at the Scale of the CA Grid
A single regular octahedron is defined as a constituent element of a grain. Each cell corresponds to a regular octahedron. The sizes of regular octahedrons are characterized by the distance (SL) between any two adjacent vertices. Its three orthogonal axes ([100]/[010]/[001]) represent the directions of dendrite tips. Meanwhile, in order to avoid the infinite growth problem of the single regular octahedron, the distance (SL) is limited to 2·△x (△x is the cell size). The grain envelope, which consists of a number of regular octahedrons, is also shown as a larger regular octahedron in the uniform temperature field. A single grain whose crystallographic orientation is assumed to be the Euler angles (φ, θ, ψ) with respect to the axis of CA, as shown in Figure 1. In this present work, the nucleated cell is defined as a solid cell. And all its neighbor cells are defined as interface cells with a solid fraction of ε (ε→0). When the solid fraction of the interface cells reaches 1 or greater than 1, all the first-order, second-order, and third-order neighbors' liquid cells of the former interface cell transform to interface cells. At this point, it is guaranteed that at least one of the 26 neighbor cells of the interface cells evolves to a solid cell. The Schematic diagram of cell state is shown in Figure 2. For the octahedron corresponding to the nucleated cell, the center of octahedron coincides with the center of nucleated cell. While the center of octahedron corresponding to the interface cell also coincides with the center of nucleated cell. Similar to 3D decentered square algorithm, the tips length of octahedron (L) from its center to any of its six vertices are given by: In this present work, the nucleated cell is defined as a solid cell. And all its neighbor cells are defined as interface cells with a solid fraction of ε (ε→0). When the solid fraction of the interface cells reaches 1 or greater than 1, all the first-order, second-order, and third-order neighbors' liquid cells of the former interface cell transform to interface cells. At this point, it is guaranteed that at least one of the 26 neighbor cells of the interface cells evolves to a solid cell. The Schematic diagram of cell state is shown in Figure 2. For the octahedron corresponding to the nucleated cell, the center of octahedron coincides with the center of nucleated cell. While the center of octahedron corresponding to the interface cell also coincides with the center of nucleated cell. Similar to 3D decentered square algorithm, the tips length of octahedron (L) from its center to any of its six vertices are given by: where t 2 represents the next time, t 1 represents the current time. The initial tip length of the grain corresponding to the nucleated cell is taken as the maximum value ( √ 2·∆x). The initial tip length of the grain corresponding to the interface cell is taken as ε (ε→0).
where t2 represents the next time, t1 represents the current time. The initial tip length of the grain corresponding to the nucleated cell is taken as the maximum value (√2·△x). The initial tip length of the grain corresponding to the interface cell is taken as ε (ε→0). During the evolution of grains, due to the displacement between the center of the interface cell and the center of the growing octahedron, the growing octahedron corresponding to the interface cell will capture the cell center of the neighbors of solid cell. When the interface cell is captured, it is then switched to the solid cell. At this moment, all the remaining liquid cells being neighbors of this new solid cell will be automatically transferred from liquid to interface cells. The new octahedron corresponding to this new interface cell is initialized. Once the tip length of the growing octahedron is equal to or greater than the maximum value (√2·△x), the growing octahedron will cease to grow. The corresponding interface cell will be switched to the solid cell regardless of the solid fraction of the interface cell.
In order to realize the self-consistency between the tip length of the growing octahedron and the solid fraction of interface cell, a quantitative method for calculating the two is described below.
Figures 3 and 4 sketch the capture rule for solid cell in a time step calculation. When the center of growing octahedron corresponding to the interface cell (for example, IJ) coincides with the center of nucleated cell, the fraction increment of the interface cell IJ is calculated by: where △L represent the increment of the tip length of the growing octahedron, L1 represent the final tip length (OA1) of the growing octahedron. During the evolution of grains, due to the displacement between the center of the interface cell and the center of the growing octahedron, the growing octahedron corresponding to the interface cell will capture the cell center of the neighbors of solid cell. When the interface cell is captured, it is then switched to the solid cell. At this moment, all the remaining liquid cells being neighbors of this new solid cell will be automatically transferred from liquid to interface cells. The new octahedron corresponding to this new interface cell is initialized. Once the tip length of the growing octahedron is equal to or greater than the maximum value ( √ 2·∆x), the growing octahedron will cease to grow. The corresponding interface cell will be switched to the solid cell regardless of the solid fraction of the interface cell.
In order to realize the self-consistency between the tip length of the growing octahedron and the solid fraction of interface cell, a quantitative method for calculating the two is described below.
Figures 3 and 4 sketch the capture rule for solid cell in a time step calculation. When the center of growing octahedron corresponding to the interface cell (for example, I J ) coincides with the center of nucleated cell, the fraction increment of the interface cell I J is calculated by: where ∆L represent the increment of the tip length of the growing octahedron, L 1 represent the final tip length (OA 1 ) of the growing octahedron. Subsequently, as soon as the center of the interface cell IJ is captured by the growing octahedron, this interface cell is immediately switched to solid cell and the remaining liquid cells adjacent to this new solid cell will be transformed into new interface cells (for example, IJK), as shown in Figure 4. As mentioned above, the new growing octahedron corresponding to the new interface Subsequently, as soon as the center of the interface cell I J is captured by the growing octahedron, this interface cell is immediately switched to solid cell and the remaining liquid cells adjacent to this new solid cell will be transformed into new interface cells (for example, I JK ), as shown in Figure 4. As mentioned above, the new growing octahedron corresponding to the new interface cell I JK is initialized. Especially, the center (O 0 ) of the new growing octahedron is also located at the center of the octahedron corresponding to the original interface cell I J . And the tip length of the new growing octahedron is initialized with the maximum value L m (L m = √ 2·∆x − ε). As a result, before the solid fraction of the new interface cell I JK reaches 1, the center (O 2 ) of the new growing octahedron has moved gradually along the direction of growth so as to maintain its continuous growth. In this case, the fraction increment of the new interface cell I JK is calculated by: where L 2 represent the final tip length (O 0 A 3 ) of the growing octahedron. Subsequently, as soon as the center of the interface cell IJ is captured by the growing octahedron, this interface cell is immediately switched to solid cell and the remaining liquid cells adjacent to this new solid cell will be transformed into new interface cells (for example, IJK), as shown in Figure 4. As mentioned above, the new growing octahedron corresponding to the new interface cell IJK is initialized. Especially, the center (O0) of the new growing octahedron is also located at the center of the octahedron corresponding to the original interface cell IJ. And the tip length of the new growing octahedron is initialized with the maximum value Lm (Lm = √2·△x − ε). As a result, before the solid fraction of the new interface cell IJK reaches 1, the center (O2) of the new growing octahedron has moved gradually along the direction of growth so as to maintain its continuous growth. In this case, the fraction increment of the new interface cell IJK is calculated by: where L2 represent the final tip length (O0A3) of the growing octahedron.

Comparison with the Benchmark Experimental Data
Firstly, to validate the numerical calculation program for mesoscopic grain growth in this study, the dendritic tip growth rate of a succinonitrile-acetone mixtures are calculated and compared with benchmark experimental data [34]. According to the microsolvability theory, the stability parameter remains unchanged within a certain range of Peclet number. Here, two stability parameters are considered. Figure 5 presents the dendritic tip growth rate as a function of initial concentration in succinonitrile-acetone mixtures. Firstly, to validate the numerical calculation program for mesoscopic grain growth in this study, the dendritic tip growth rate of a succinonitrile-acetone mixtures are calculated and compared with benchmark experimental data [34]. According to the microsolvability theory, the stability parameter remains unchanged within a certain range of Peclet number. Here, two stability parameters are considered. Figure 5 presents the dendritic tip growth rate as a function of initial concentration in succinonitrile-acetone mixtures. As shown in Figure 5, under both undercooling conditions of 0.5 K and 0.9 K, the dendrite tip growth rate calculated by numerical methods increases first and then decreases with the initial concentration, which is consistent with the experimental results. For the initial concentration range of 0-0.25 (mol. %), the curve of the dendritic tip growth rate with a stability parameter of 0.025 is approximately matched with the experimental results. For the initial concentration range of 0.25-0.45 (mol%), the curve of the dendritic tip growth rate with a stability parameter of 0.0195 is approximately matched with the experimental results. It is confirmed that the stability parameters match the experimental results only in a certain range of Peclet numbers.

Comparison with the Analytical Model of Grain Growth
The model and algorithm validation for grain growth simulation is made by comparing the numerical results with the analytical model predictions [14,36]. In the analytical growth model, it was assumed that the grain growth was controlled by the solute diffusion and only applicable to the low Peclet number regime. This is basically similar to the assumption of the LGK model. There are two kinds of grain branches predicted by the analytical model, here the external envelope is assumed to be the grain envelope. Figure 6 presents the schematic diagram of two-dimensional dendrite growth envelope trajectory under directional temperature. In the process of two-dimensional dendritic growth under directional solidification, the coordinates of the external envelope trajectory can be calculated by Equation (11) [36]: As shown in Figure 5, under both undercooling conditions of 0.5 K and 0.9 K, the dendrite tip growth rate calculated by numerical methods increases first and then decreases with the initial concentration, which is consistent with the experimental results. For the initial concentration range of 0-0.25 (mol. %), the curve of the dendritic tip growth rate with a stability parameter of 0.025 is approximately matched with the experimental results. For the initial concentration range of 0.25-0.45 (mol%), the curve of the dendritic tip growth rate with a stability parameter of 0.0195 is approximately matched with the experimental results. It is confirmed that the stability parameters match the experimental results only in a certain range of Peclet numbers.

Comparison with the Analytical Model of Grain Growth
The model and algorithm validation for grain growth simulation is made by comparing the numerical results with the analytical model predictions [14,36]. In the analytical growth model, it was assumed that the grain growth was controlled by the solute diffusion and only applicable to the low Peclet number regime. This is basically similar to the assumption of the LGK model. There are two kinds of grain branches predicted by the analytical model, here the external envelope is assumed to be the grain envelope. Figure 6 presents the schematic diagram of two-dimensional dendrite growth envelope trajectory under directional temperature. In the process of two-dimensional dendritic growth under directional solidification, the coordinates of the external envelope trajectory can be calculated by Equation (11) [36]: where V L is the moving rate of liquid isotherm, ∆T L is the dendrite tip undercooling, ∆T 0 is the undercooling at initial time t 0 , θ is the dendrite orientation angle, A is a constant coefficient in LGK model, f (θ), g(θ) and S <10> are directional functions related to the orientation angle respectively, which are given by Equations (12)- (14).
where VL is the moving rate of liquid isotherm, △TL is the dendrite tip undercooling, △T0 is the undercooling at initial time t0, θ is the dendrite orientation angle, A is a constant coefficient in LGK model, f(θ), g(θ) and S<10> are directional functions related to the orientation angle respectively, which are given by Equations (12)- (14).
Since the above analytical model is based on two-dimensional dendrite growth, it is necessary to extend this method to three-dimensional grain growth. In this present study, the projection method and the three cosine theorem are used to calculate the outermost envelope trajectory of a three-dimensional grain. Besides, since the analytical model was restricted to a quadratic approximation of the velocity-undercooling relationship, the growth rate was temporarily for the [ Since the above analytical model is based on two-dimensional dendrite growth, it is necessary to extend this method to three-dimensional grain growth. In this present study, the projection method and the three cosine theorem are used to calculate the outermost envelope trajectory of a three-dimensional grain. Besides, since the analytical model was restricted to a quadratic approximation of the velocity-undercooling relationship, the growth rate was temporarily simplified as the same quadratic approximation relationship of local undercooling in the numerical simulation. Figure 7 presents the 3D envelopes and tips length of a single grain growing in a uniform temperature field and in a Bridgman-like situation. In Figure 7a, the cooling rate is set to 0 K/s and the thermal gradient is set to 0 K/m, the tips length (L 1 , L 2 , L 3 , L 4 , L 5 , L 6 ) of the grain corresponding to Figure 7a are shown in Figure 7b. In Figure 7c, the cooling rate is set to −0.1 K/s and the thermal gradient along the Z-axis is set to 250 K/m, the tips length (L 1 , L 2 , L 3 , L 4 , L 5 , L 6 ) of the grain corresponding to Figure 7c are shown in Figure 7d. The cell size ∆x is set to 0.0001 m.
simulation. Figure 7 presents the 3D envelopes and tips length of a single grain growing in a uniform temperature field and in a Bridgman-like situation. In Figure 7a, the cooling rate is set to 0 K/s and the thermal gradient is set to 0 K/m, the tips length (L1, L2, L3, L4, L5, L6) of the grain corresponding to Figure 7a are shown in Figure 7b. In Figure 7c, the cooling rate is set to −0.1 K/s and the thermal gradient along the Z-axis is set to 250 K/m, the tips length (L1, L2, L3, L4, L5, L6) of the grain corresponding to Figure 7c are shown in Figure 7d. The cell size △x is set to 0.0001 m. As shown in Figure 7a,c the 3D views of analytical and numerical predictions, for a given crystallographic orientation, are approximately regular octahedrons in a uniform temperature field. In Bridgman-like situations, due to the effects of deep undercooling, the tips length (L4, L6) of the grain closest to the negative thermal gradient direction are the longer distances. Altogether, the 3D grain envelopes (a2, c2) simulated by the developed model and algorithm in both cases are visually reproduced and coincident with the shape (a1, c1) predicted by the analytical model. Furthermore, from the Figures 7b and 6d, it can be noted that there are only a few mesh size gaps among the tips length of the grain in both predictions, which also approximately present the same distribution. As shown in Figure 7a,c the 3D views of analytical and numerical predictions, for a given crystallographic orientation, are approximately regular octahedrons in a uniform temperature field. In Bridgman-like situations, due to the effects of deep undercooling, the tips length (L 4 , L 6 ) of the grain closest to the negative thermal gradient direction are the longer distances. Altogether, the 3D grain envelopes (a2, c2) simulated by the developed model and algorithm in both cases are visually reproduced and coincident with the shape (a1, c1) predicted by the analytical model. Furthermore, from the Figure 7b,d, it can be noted that there are only a few mesh size gaps among the tips length of the grain in both predictions, which also approximately present the same distribution.

Model Application to Nickel-Based Superalloys
The present model is used to simulate the formation of equiaxed and columnar dendritic grains during nickel-based superalloys solidification process. The elemental composition (mass fraction, wt.%) of the material is Ni-7.82Cr-5.34Co-2.25Mo-4.88W-6.02Al-1.94Ti-3.49Ta. For multicomponent alloys, a large number of simplifications and calculations are required due to the complex multicomponent coupling solidification process. Here, a multicomponent pseudo-binary alloy method [37] is used. The material properties and model parameters used in these simulations are listed in Table 1. 3.0 × 10 −9 Solute diffusion coefficient in the solid (m 2 ·s −1 ) 3.0 × 10 −12

Equiaxed Dendritic Grains
The calculation domain is uniformly divided into 200 × 200 × 200 meshes with mesh size of 0.0001 m. It is assumed that the temperature field is homogeneous and cooled down from liquidus temperature at different cooling rates leading to nucleation and growth process. The cooling rates are set to 1 K·s −1 , 10 K·s −1 , and 100 K·s −1 , respectively. Here, the nucleation process only takes into account the nucleation inside the melt. It must be stressed that the nucleation parameters are affected by many factors, such as random nucleation, solidification conditions, etc. It cannot be uniquely determined or calculated in advance. For the ease of numerical calculations, the nucleation parameter is assumed to be known in advance. In the bulk of the melt, the total density of grains is 1 × 10 10 m −3 , the median undercooling is 3 K, and the standard deviation of undercooling is 0.5 K. Figure 8 presents the size distribution of equiaxed dendritic grains at different cooling rates.
The present model is used to simulate the formation of equiaxed and columnar dendritic grains during nickel-based superalloys solidification process. The elemental composition (mass fraction, wt%) of the material is Ni-7.82Cr-5.34Co-2.25Mo-4.88W-6.02Al-1.94Ti-3.49Ta. For multicomponent alloys, a large number of simplifications and calculations are required due to the complex multicomponent coupling solidification process. Here, a multicomponent pseudo-binary alloy method [37] is used. The material properties and model parameters used in these simulations are listed in Table 1. 1.0 × 10 −7 Solute diffusion coefficient in the liquid (m 2 ·s −1 ) 3.0 × 10 −9 Solute diffusion coefficient in the solid (m 2 ·s −1 ) 3.0 × 10 −12

Equiaxed Dendritic Grains
The calculation domain is uniformly divided into 200 × 200 × 200 meshes with mesh size of 0.0001 m. It is assumed that the temperature field is homogeneous and cooled down from liquidus temperature at different cooling rates leading to nucleation and growth process. The cooling rates are set to 1 K·s −1 , 10 K·s −1 , and 100 K·s −1 , respectively. Here, the nucleation process only takes into account the nucleation inside the melt. It must be stressed that the nucleation parameters are affected by many factors, such as random nucleation, solidification conditions, etc. It cannot be uniquely determined or calculated in advance. For the ease of numerical calculations, the nucleation parameter is assumed to be known in advance. In the bulk of the melt, the total density of grains is 1 × 10 10 m −3 , the median undercooling is 3 K, and the standard deviation of undercooling is 0.5 K. Figure 8 presents the size distribution of equiaxed dendritic grains at different cooling rates.  Figure 8 shows that, with the increase of cooling rate, the grain size of equiaxed grains decreases obviously, and the grain density increases from 0.13 mm −2 to 4.08 mm −2 . From the data analysis of the two-dimensional section, it can be seen that the number of grains approximately  Figure 8 shows that, with the increase of cooling rate, the grain size of equiaxed grains decreases obviously, and the grain density increases from 0.13 mm −2 to 4.08 mm −2 . From the data analysis of the two-dimensional section, it can be seen that the number of grains approximately presents a normal distribution curve at d/<d> from 0 to 2.0. With higher cooling rates the maximum of the grain number as function of d/<d> grows, while its dispersion in grain-size lowers. The highest grain number-density is always given close to the median grain-size. It is indicated that the increase of cooling rate can significantly promote the grain refinement. These results are consistent with the experimental results of solidification microstructures [38], which preliminarily validate the effectiveness of the proposed model and algorithm.

Columnar Dendritic Grains
Cooling rate and temperature gradient are two important parameters for directional solidification microstructure. Since the thermal diffusion rate is more than several orders of magnitude larger than the solute diffusion rate, it can be assumed that the temperature is uniform at the scale of a grain. At a fixed temperature gradient (G) and cooling rate (C) the temperature isotherms of the whole calculation domain will change at the same drawing speed (R = C/G). The calculation domain is uniformly divided into 60 × 180 × 600 cubic cells with size of 0.0001 m. Two types of nucleation parameters are considered. At the bottom surface, the total density of grains is 5.0 × 10 5 m −3 , the mean undercooling is 0 K, and the standard deviation is 0 K. In the bulk of the melt, the total density of grains is 1 × 10 10 m −3 , the mean undercooling is 3 K, and the standard deviation is 0.5 K.

Equiaxed-to-Columnar Transformation (ECT) for Different Cooling Rates
The temperature gradient is set to 1000 K·m −1 . And the cooling rates are set to 1 K·s −1 , 10 K·s −1 , and 100 K·s −1 , respectively. Figure 9 presents the grain size distributions for different cooling rates. presents a normal distribution curve at d/<d> from 0 to 2.0. With higher cooling rates the maximum of the grain number as function of d/<d> grows, while its dispersion in grain-size lowers. The highest grain number-density is always given close to the median grain-size. It is indicated that the increase of cooling rate can significantly promote the grain refinement. These results are consistent with the experimental results of solidification microstructures [38], which preliminarily validate the effectiveness of the proposed model and algorithm.

Columnar Dendritic Grains
Cooling rate and temperature gradient are two important parameters for directional solidification microstructure. Since the thermal diffusion rate is more than several orders of magnitude larger than the solute diffusion rate, it can be assumed that the temperature is uniform at the scale of a grain. At a fixed temperature gradient (G) and cooling rate (C) the temperature isotherms of the whole calculation domain will change at the same drawing speed (R = C/G). The calculation domain is uniformly divided into 60 × 180 × 600 cubic cells with size of 0.0001 m. Two types of nucleation parameters are considered. At the bottom surface, the total density of grains is 5.0 × 10 5 m −3 , the mean undercooling is 0 K, and the standard deviation is 0 K. In the bulk of the melt, the total density of grains is 1 × 10 10 m −3 , the mean undercooling is 3 K, and the standard deviation is 0.5 K.

Equiaxed-to-Columnar Transformation (ECT) for Different Cooling Rates
The temperature gradient is set to 1000 K·m −1 . And the cooling rates are set to 1 K·s −1 , 10 K·s −1 , and 100 K·s −1 , respectively. Figure 9 presents the grain size distributions for different cooling rates.  Figure 9 shows, that for different cooling rates, the developed solidification microstructures of the calculated domain gradually changes (with respect to the size distribution of grains). At the cooling rate of 1 K·s −1 , some thick columnar dendritic grains are shown in the calculation domain, and the number of grains decreases gradually along the Z axis. At the cooling rate of 10 K·s −1 , the grain structure in the calculation domain is composed of both columnar and fine equiaxed dendritic grains, and the number of grains increases first and then decreases. At the cooling rate of 100 K·s −1 , almost all grains in the calculation domain are small equiaxed dendritic grains, and the number of grains increases suddenly and decreases slowly with increasing Z. From the above results, it can be inferred that at low cooling rate the grain growth rate of directional solidification is faster than the grain nucleation rate in the supercooled region of melt. In addition, due to the boundary limitation of the calculation domain, some grains deviating from the direction of temperature gradient will be eliminated at the boundary, resulting in the gradual decrease of the number of grains along the Z axis.

Equiaxed-to-Columnar Transformation (ECT) for Different Temperature Gradients
The temperature gradient is set to 500 K·m −1 , 1000 K·m −1 , and 1500 K·m −1 , respectively. The cooling rates are set to 10 K·s −1 . Figure 10 presents the grain size distributions for different temperature gradients.  Figure 9 shows, that for different cooling rates, the developed solidification microstructures of the calculated domain gradually changes (with respect to the size distribution of grains). At the cooling rate of 1 K·s −1 , some thick columnar dendritic grains are shown in the calculation domain, and the number of grains decreases gradually along the Z axis. At the cooling rate of 10 K·s −1 , the grain structure in the calculation domain is composed of both columnar and fine equiaxed dendritic grains, and the number of grains increases first and then decreases. At the cooling rate of 100 K·s −1 , almost all grains in the calculation domain are small equiaxed dendritic grains, and the number of grains increases suddenly and decreases slowly with increasing Z. From the above results, it can be inferred that at low cooling rate the grain growth rate of directional solidification is faster than the grain nucleation rate in the supercooled region of melt. In addition, due to the boundary limitation of the calculation domain, some grains deviating from the direction of temperature gradient will be eliminated at the boundary, resulting in the gradual decrease of the number of grains along the Z axis.

Equiaxed-to-Columnar Transformation (ECT) for Different Temperature Gradients
The temperature gradient is set to 500 K·m −1 , 1000 K·m −1 , and 1500 K·m −1 , respectively. The cooling rates are set to 10 K·s −1 . Figure 10 Figure 10 shows for the calculated domain, that with increasing temperature gradient, the solidification microstructure develops to the opposite of the results presented in Figure 8. In other words, at the temperature gradient of 500 K·m −1 , almost all grains in the calculation domain are small equiaxed dendritic grains, and the number of grains increases suddenly and decreases slowly. At the temperature gradient of 1500 K·m −1 , some thick columnar dendritic grains are shown in the calculation domain, and the number of grains decreases gradually along the Z axis. It can also be deduced that a rapid growth of columnar dendritic grains is promoted more by an increasing temperature gradient than by an increasing nucleation rate.

Columnar-to-Equiaxed Transformation (CET)
The ECT and CET are common phenomena of the grain evolution in solidification process. Under different solidification conditions, the transformations rule of the two influences the grain size distribution in the casting. During directional solidification, the ECT transformation is beneficial to the growth of columnar dendritic grains, and it is also affected by the nucleation parameters and cooling rate. Since the cooling rate and the temperature gradient are difficult to control, the casting  Figure 10 shows for the calculated domain, that with increasing temperature gradient, the solidification microstructure develops to the opposite of the results presented in Figure 8. In other words, at the temperature gradient of 500 K·m −1 , almost all grains in the calculation domain are small equiaxed dendritic grains, and the number of grains increases suddenly and decreases slowly. At the temperature gradient of 1500 K·m −1 , some thick columnar dendritic grains are shown in the calculation domain, and the number of grains decreases gradually along the Z axis. It can also be deduced that a rapid growth of columnar dendritic grains is promoted more by an increasing temperature gradient than by an increasing nucleation rate.

Columnar-to-Equiaxed Transformation (CET)
The ECT and CET are common phenomena of the grain evolution in solidification process. Under different solidification conditions, the transformations rule of the two influences the grain size distribution in the casting. During directional solidification, the ECT transformation is beneficial to the growth of columnar dendritic grains, and it is also affected by the nucleation parameters and cooling rate. Since the cooling rate and the temperature gradient are difficult to control, the casting of the metallic alloy may exhibit a non-uniform distribution of the temperature field. Especially, when the temperature gradient gradually decreases to a negative temperature gradient, this will lead to the CET transformation. To characterize the positive-to-negative transition in the temperature gradient, the temperature gradient is set to 1000 K·m −1 from mesh point 0 to 300 in the direction of the Z axis, and the temperature gradient is set to −300 K·m −1 from mesh point 300 to 600 in the direction of the Z axis. The cooling rate is set to 10 K·s −1 . The other grain nucleation and growth parameters set here are the same as those above Section 3.3.2. Figure 11 presents the simulation results of the CET transformation. of the metallic alloy may exhibit a non-uniform distribution of the temperature field. Especially, when the temperature gradient gradually decreases to a negative temperature gradient, this will lead to the CET transformation. To characterize the positive-to-negative transition in the temperature gradient, the temperature gradient is set to 1000 K·m −1 from mesh point 0 to 300 in the direction of the Z axis, and the temperature gradient is set to −300 K·m −1 from mesh point 300 to 600 in the direction of the Z axis. The cooling rate is set to 10 K·s −1 . The other grain nucleation and growth parameters set here are the same as those above Section 3.2.2. Figure 11 presents the simulation results of the CET transformation. As shown in Figure 11, in the calculation domain, the grains first nucleate and grow at the bottom surface. From mesh point 0 to 300 in the direction of the Z axis, under the positive temperature gradient, grains begin to compete with each other along the Z axis, which is similar to the ECT transformation process mentioned above. From mesh point 300 to 600 in the direction of the Z axis, the grains nucleate and grow at the top surface under the negative temperature gradient. Compared with the ECT transformation process, the grains grow in the form of a large number of equiaxed crystals. Finally, the two parts of the grain stop growing at intermediate position, thus forming the CET transformation process. The above results show a good agreement with the results As shown in Figure 11, in the calculation domain, the grains first nucleate and grow at the bottom surface. From mesh point 0 to 300 in the direction of the Z axis, under the positive temperature gradient, grains begin to compete with each other along the Z axis, which is similar to the ECT transformation process mentioned above. From mesh point 300 to 600 in the direction of the Z axis, the grains nucleate and grow at the top surface under the negative temperature gradient. Compared with the ECT transformation process, the grains grow in the form of a large number of equiaxed