Numerical Simulation of Three-Dimensional Dendrite Movement Based on the CA–LBM Method

: At present, the calculation of three-dimensional (3D) dendrite motion using the cellular automata (CA) method is still in its infancy. In this paper, a 3D dendrite motion model is constructed. The heat, mass, and momentum transfer process in the solidiﬁcation process of the alloy melt are calculated using a 3D Lattice–Boltzmann method (LBM). The growth process for the alloy microstructure is calculated using the CA method. The interactions between dendrites and the melt are assessed using the Ladd method. The solid–liquid boundary of the solute ﬁeld in the movement process is assessed using the solute extrapolation method. The translational velocity of the equiaxed crystals in motion is calculated using the classical mechanical law. The rationality of the model is veriﬁed and the movement of single and multiple 3D equiaxed crystals is simulated. Additionally, the difference between 3D dendrite movement and two-dimensional (2D) dendrite movement is analyzed. The results demonstrate that the growth of moving dendrites is asymmetric. The growth velocity and falling velocity of the dendrite in the 3D model are faster than that in 2D model, while the simulation result is more realistic than that of the 2D model. When multiple dendrites move, the movement direction of the dendrites will change due to the merging of ﬂow ﬁelds and other factors.


Introduction
In recent years, with the increase in research on microstructure simulation and the rapid development of computer technology, scholars have realized that the movement of dendrites will change their growth mode, and have begun to study the influence of dendrite movement more deeply. As early as 1998, Appolaire et al. [1] studied the influence of the settlement process of equiaxed grains on dendrite growth. They facilitated NH 4 Cl-H 2 O dendrite growth in a tube filled with supercooled solution. The dendrites began to settle down under the action of gravity. They then photographed the settling process for NH 4 Cl equiaxed crystal and analyzed the changes in dendrite size and settling velocity with time. However, due to the limitations of experimental methods and the rapid development of computer technology, scholars have begun to use numerical simulation methods.
In 2007, Amerg et al. [2] used the phase field method (PF) model to simulate the melt flow and convective heat transfer process, as well as the growth and fall of dendrites under the action of gravity, revealing the complex coupling relationship between the settlement and evolution of grains. Badillo et al. [3] observed and measured the settlement of equiaxed grains in alloy succinonitrile-acetone and established a prediction model for solid fraction ( f s ) evolution under convection conditions based on thermal and solute equilibriums. Wu et al. [4] proposed a modified volume average equiaxed solidification model and verified it. In this model, nucleation, dendrite growth, grain movement, and melt convection are considered, and the effects of these factors on the macrosegregation of the dendrite morphology are analyzed. In 2012, Karagadde et al. [5] established a new model to simulate the growth and motion of equiaxed crystals. By combining the volume of fluid (VOF) method, immersed boundary method (IBM), and enthalpy method, the force and growth behavior of dendrites under the action of the flow field were studied. The results demonstrated that the flow field can change the position and growth mode of dendrites. In 2013, Medvedev [6] proposed a method to simulate the solidification of dendrites driven by solutes. The liquid flow and motion of single equiaxed crystals in the melt were studied using the LBM coupled with the multi-PF method, demonstrating that there may be movement and rotation of the solid phase in the melt. Kharicha et al. [7] studied the solution flow and columnar solidification process in NH 4 Cl-H 2 O ingots using dual particle image velocimetry (PIV). The kinetic energy of flow was calculated using PIV, while the interaction between equiaxed crystals and the melt flow was analyzed according to the relative velocity. It was found that the settling velocity of a single equiaxed crystal is 41 times smaller than that of a spherical crystal of the same size. It was demonstrated that the complex shape of equiaxed crystals affects the solid-liquid interaction, and then affects the settling velocity. Rojas et al. [8] used phase-field LBM to simulate the growth and movement of single equiaxed crystals in alloy melt under shear flow. In the same year, Takaki and Rojas et al. [9] improved the calculation efficiency by introducing GPU calculation technology and used this model to simulate the growth and movement of equiaxed grains in a binary alloy single-phase solid solution. Ludwig et al. [10] described the settlement and deposition of crystals in the melt using the combined Eulerian-Eulerian volume-averaged approach. In 2017, Xinbo et al. [11] calculated the settlement of single dendrites and multidendrites in the melt of an aluminum-based alloy. The results demonstrated that the moving and stationary dendrites, considering natural convection, form completely different morphologies. In 2018, Takaki et al. [12] established the LBM coupled PF model to simulate the movement and collision of multiple grains, then used the model to simulate the morphological evolution process for multiple grains under the action of a complex flow field; however, due to the complexity of the calculation process, the PF method can only calculate a small number of equiaxed grains, so the current PF method is not suitable for calculating the solidification process for ingots containing a large number of grains [13].
The CA method can efficiently calculate the growth process for dendrites in large-scale computational domains. The calculation accuracy is close to that of the PF method, as indicated by the calculation results in [14]. There have been many studies on the calculation of microstructure morphologies using the CA method; however, the calculation of the movement behavior of equiaxed crystals using the CA method is still in its infancy. At present, only Liu et al. [15] and Bai et al. [16] have used the 2D CA and LBM methods to simulate the movement of grains in the melt and the vertical drop of multiple grains. However, it is difficult for the 2D model to reflect the 3D dendrite's structure and movement process. Therefore, a CA-LBM coupling model is proposed in this paper to simulate the settlement of 3D equiaxed crystals. The CA method is used to calculate the growth of dendrites, while the LBM is used to establish the macrofield transport model. The Ladd method [17] is used to calculate the solid-liquid boundary and the solute extrapolation distribution method is used to distribute the solute at the solid-liquid boundary, while Newton's law of mechanics is used to calculate the settling velocity in the process of dendrite movement, so as to calculate the 3D dendrite movement.

Materials and Methods
In this study, Al-4.7% Cu alloy was selected as the research material, with the physical property parameters presented in Table 1.

LBM Model
In this study, the D 3 Q 15 model was used to calculate the temperature field, flow field, and solute field in the process of 3D dendrite movement and growth, while the D 2 Q 9 model was used to calculate the three fields in the process of 2D dendrite movement and growth [18].
The distribution functions of the flow field in LBM are as follows: where e i is the particle velocity in direction i, ∆t is the time step, and τ f is the flow field relaxation time. Here, f i (x, t) is the distribution function, which represents the probability of particles appearing in position x, while time t. f eq i (x, t) is the equilibrium distribution function, which can be calculated by: where c is the lattice velocity and ω i is the weight coefficient; according to different models, ω i takes different values.
The relationship between τ f and viscosity ν can be expressed as: The resultant force F of the particle can be expressed as: where ρ 0 is the initial density of the fluid, T is the current temperature, T 0 is the initial temperature, C is the current concentration, C 0 is the initial concentration, β T is the volume expansion coefficient of the temperature change, β C is the volume expansion coefficient of the concentration change, and g is the gravity acceleration. The density and velocity of the macroscopic physical quantity can be evaluated by: where M = 14 in D3Q15 and M = 8 in D2Q9. The distribution functions of the temperature field and solute field are as follows: The equilibrium distribution functions of the temperature field and solute field are: In Equations (7) and (8), H i and G i respectively represent the latent heat and concentration source terms at the solid-liquid interface due to latent heat release and solute redistribution: where L is the latent heat of solidification, c p is the specific heat capacity, C l is the composition in the liquid phase, k is the equilibrium distribution coefficient, and ∆ f s is the increase in the solid fraction of an interface cell in a time step. The relationships between the temperature field relaxation time τ α and the thermal diffusion coefficient α, concentration field relaxation time τ D , and solute diffusion coefficient D are as follows: The macroscopic temperature and concentration of the fluid can be obtained by summing the corresponding distribution functions: In the process of solidification, the solute will be discharged from the growing cells to the surrounding liquid cells. In order to conserve the solute, the modified solute distribution model proposed by Li [19] is used in the calculation of the solute redistribution. In this model, it is assumed that the diffusion only occurs at the liquid phase and solidliquid interface. It is assumed that the growth cell is named x, the initial solid concentration is C 0 s , the solid fraction is f 0 s , and the solute concentration in the liquid phase is C 0 l ; then, it is judged whether there are liquid cells around the growing cell x. When there are liquid phase cells around cell x and the solid fraction increment ∆ f s ≥ 1 − f 0 s , the liquid and solid phase concentration of cell x are: The amount of solute discharged during cell x growth can be expressed as: When there are liquid phase cells around cell x and ∆ f s < 1 − f 0 s ,the liquid phase concentration of cell x is: If C l > C eq l , let the value of C l be equal to C eq l . The amount of solute discharged during cell x growth is: However, if C l ≤ C eq l , then ∆C = 0, and the solid phase concentration of cell x is: When there is no liquid phase cell in the surrounding cell of cell x and ∆ f s ≥ 1 − f 0 s , the solute discharged amount ∆C = 0, the liquid phase concentration in cell x C l = 0 , and the solid phase concentration are: When there is no liquid phase cell in the surrounding cells of cell x, and ∆ f s < 1 − f 0 s , the solute discharged amount ∆C = 0 and the liquid and solid phase concentrations are respectively: If there are liquid phase cells around cell x, the discharged solute is weighted into the liquid phase cells.

CA Model
In this study, the solute diffusion model proposed by Zhu [20] was used to calculate the dendrite growth. The solid phase increment of the growth cell in this model can be calculated by: where k is the equilibrium partition coefficient, C l is the actual solute concentration of the interfacial cell, and C eq l is the equilibrium solute concentration of interfacial cell. The calculation details are provided in [20].
In order to eliminate the grid-induced anisotropy caused by cell capture as much as possible, we selected the 26 neighborhood capture mode to calculate 3D dendrites and the 8 neighborhood capture mode to calculate 2D dendrites; that is, in the 3D calculation, the interface cells capture 26 surrounding liquid cells, while in the 2D calculation, the interface cells capture 8 surrounding liquid cells.

Boundary Treatment
The boundary conditions in this paper include computational domain boundary conditions and solid-liquid interface boundary conditions. The rebound scheme is used at the outer boundary of the flow field, while the half step rebound scheme is used at the solidliquid boundary [17]. For the temperature field, the non-equilibrium extrapolation scheme is adopted for the boundary of the calculation domain. Because the thermal diffusion coefficients of solid and liquid phases are equal and constant, no special treatment is used for the solid-liquid boundary. For the solute field, the boundary of the calculation domain adopts the rebound scheme; however, the Ladd method cannot be used because of the large difference of the concentration diffusion coefficients between the solid and liquid phases. In this paper, the solute extrapolation distribution method [16] is used to deal with the solid-liquid boundary.

Verification
The rationality test used for calculation of the CA model was confirmed in relevant studies [21,22], while the calculation of dendrite growth using the LBM-CA model was verified by Pian [23], which will not be described in this paper. The calculation of the interaction force at the solid-liquid interface and the degree of artificial error introduced in the movement of the dendrite need to be verified, which will be verified by comparing the analytical solution of the settling velocity of the 3D spherical particles in an infinite pipe with the calculation results and by calculating the dendrite changes after translation.

Settlement of Small, Spherical 3D Particles in Infinite Length Pipeline
Glowiski [24] proposed that when the physical parameters of a ball are fixed, the final stable settlement velocity of the ball in the infinite pipe is also fixed. When the fluid is at a low Reynolds number, the influence of viscous force on the flow field is greater than that of inertial force and the disturbance of the flow velocity in the flow field will be attenuated due to the viscous force, meaning the fluid flow will be relatively stable.
The final settling velocity of spherical particles can be calculated by [15]: where D is the diameter of the spherical particle, ρ s and ρ l are the densities of the spherical particles and fluid, respectively, g is the acceleration of gravity, and η is the dynamic viscosity of the fluid. Here, K is the correction factor of the resistance of the pipe wall to the particle, representing the influence of the long pipe wall on the particle: In the simulation, the calculation domain is divided into 80 × 80 × 500 cubic grids, the grid size h is 0.003 cm, the diameter of the 3D spherical particle is 0.048 cm, and the width between plates L is 0.24 cm. The density of the fluid is ρ l = 1.0 g/cm 3 , the density of solid particle is ρ s = 1.02 g/cm 3 , the acceleration of gravity is g = 980.0 cm/s 2 , and the viscosity of the fluid is η = 0.33 g/(cm·s). At the initial time, the static spherical particle is placed at the location indicated in Figure 1. The velocity at the entrance boundary is 0, and the velocity derivative at the exit boundary is 0. The non-equilibrium extrapolation scheme is used to deal with six wall boundaries. When the calculation time t  = 3.125 × 10 −5 s, the maximum Reynolds number is 469.37, which is in good agreement with the results from the studies by Wan [25], Uhlmann [26], and Glowiski [24]. In this paper, the spherical particle is always located on the central axis of the pipe. The spatial discretization scheme is highly symmetrical, and the numerical disturbance caused by the wall resistance is small. Figure 2 demonstrates that the spherical particle falls from the static state, then, because the gravity is greater than the resultant force of buoyancy and resistance, the particle continues to accelerate, while the resistance also increases and the acceleration decreases. Finally, the three forces of gravity, buoyancy, and resistance are balanced. The velocity is consistent with the theoretical value when it is stable, demonstrating that the treatment of When the calculation time ∆t = 3.125 × 10 −5 s, the maximum Reynolds number is 469.37, which is in good agreement with the results from the studies by Wan [25], Uhlmann [26], and Glowiski [24]. In this paper, the spherical particle is always located on the central axis of the pipe. The spatial discretization scheme is highly symmetrical, and the numerical disturbance caused by the wall resistance is small. Figure 2 demonstrates that the spherical particle falls from the static state, then, because the gravity is greater than the resultant force of buoyancy and resistance, the particle continues to accelerate, while the resistance also increases and the acceleration decreases. Finally, the three forces of gravity, buoyancy, and resistance are balanced. The velocity is consistent with the theoretical value when it is stable, demonstrating that the treatment of the moving boundary of the solid-liquid interface is reasonable. When the calculation time t  = 3.125 × 10 −5 s, the maximum Reynolds number is 469.37, which is in good agreement with the results from the studies by Wan [25], Uhlmann [26], and Glowiski [24]. In this paper, the spherical particle is always located on the central axis of the pipe. The spatial discretization scheme is highly symmetrical, and the numerical disturbance caused by the wall resistance is small. Figure 2 demonstrates that the spherical particle falls from the static state, then, because the gravity is greater than the resultant force of buoyancy and resistance, the particle continues to accelerate, while the resistance also increases and the acceleration decreases. Finally, the three forces of gravity, buoyancy, and resistance are balanced. The velocity is consistent with the theoretical value when it is stable, demonstrating that the treatment of the moving boundary of the solid-liquid interface is reasonable.

Figure 2.
Settling velocity for the spherical particle.

Verification of the Translation of A Single Dendrite
The computational domain is divided into 110 × 110 × 110 cubic grids with a size of 1 μm. The initial composition of the liquid phase 0 4.7% C  . The left boundary of the simulation domain is a constant entrance velocity 0 0.08 cm/s  u , the right boundary is a fully developed boundary, and the other four boundaries are sliding boundaries. At the initial time, a nucleus is placed at (50, 20, 50), the s f of this seed is 0.2, and the preferential growth angle is 0 rad. At the beginning, the effect of the flow field on dendrite growth is ignored, so that the dendrite grows in a static state. When the number of time steps is 1000, the growth of the dendrite is ignored, and a uniform flow field is added into the domain. The dendrite begins to move horizontally due to the flow field.
Dendrites accelerate along the horizontal direction with the pipeline flow, and the resistance of the dendrites increases as the velocity of the dendrites increases. When the drag is equal to the traction force, the dendrite will move at a uniform velocity in the horizontal direction. In order to quantitatively describe the velocity of the dendrite, the

Verification of the Translation of a Single Dendrite
The computational domain is divided into 110 × 110 × 110 cubic grids with a size of 1 µm. The initial composition of the liquid phase C 0 = 4.7%. The left boundary of the simulation domain is a constant entrance velocity u 0 = 0.08 cm/s, the right boundary is a fully developed boundary, and the other four boundaries are sliding boundaries. At the initial time, a nucleus is placed at (50, 20, 50), the f s of this seed is 0.2, and the preferential growth angle is 0 rad. At the beginning, the effect of the flow field on dendrite growth is ignored, so that the dendrite grows in a static state. When the number of time steps is 1000, the growth of the dendrite is ignored, and a uniform flow field is added into the domain. The dendrite begins to move horizontally due to the flow field.
Dendrites accelerate along the horizontal direction with the pipeline flow, and the resistance of the dendrites increases as the velocity of the dendrites increases. When the drag is equal to the traction force, the dendrite will move at a uniform velocity in the horizontal direction. In order to quantitatively describe the velocity of the dendrite, the movement velocity-time curve of dendrite is provided in Figure 3, which demonstrates that when the dendrite reaches a stable state, the velocity is the same as that of the fluid.
Crystals 2021, 11, x FOR PEER REVIEW movement velocity-time curve of dendrite is provided in Figure 3, which demo that when the dendrite reaches a stable state, the velocity is the same as that of the The change degree of the dendrite morphology is quantitatively described b ducing  : Here, 0 s f is the solid fraction before the dendrite movement and 1 s f is th fraction after the dendrite movement;  is approximately 0 at any time, which tatively demonstrates that the solid fraction will not be lost due to the translation. The change degree of the dendrite morphology is quantitatively described by introducing λ: Here, f 0 s is the solid fraction before the dendrite movement and f 1 s is the solid fraction after the dendrite movement; λ is approximately 0 at any time, which quantitatively demonstrates that the solid fraction will not be lost due to the translation.

Results and Discussions
The calculations in this paper are based on the following assumptions: (1) The melt is considered as an incompressible Newtonian fluid. (2) The temperature of the six walls in the computational domain is 913 K and the initial undercooling is 4 K. (3) The thermal conductivity of the solid phase and liquid phase is the same and the heat transfer boundary condition is adiabatic. (4) Solute diffusion only occurs in the liquid phase; there is no solute diffusion in the solid phase. The solid phase is homogeneous, and the mass transfer boundary condition is a non-diffusion boundary condition. (5) The impacts of dendrite collisions are ignored.

Comparison of Simulation Results for 3D and 2D Models under Moving Conditions
In 3D calculationg, the computational domain is divided into 140 × 140 × 140 cubic grids with a size of 0.5 µm. At the beginning of the solidification, a seed with an initial preferential growth angle of 0 rad is placed at the upper position (70, 70, 100) of the calculation domain and the nucleus grows along the coordinate axis. The initial composition of the liquid phase is C 0 = 4.7%, while k is the equilibrium distribution coefficient; therefore, the solid solute composition of the nucleus is kC 0 and the Prandtl number is 0.7. The density of the solid phase of the Al-4.7% Cu alloy is higher than that of the liquid phase, and the primary dendrites formed by solidification will grow while settling in the melt, consistent with the experimental observations [27,28]. The 2D simulation calculation domain is 400 × 800 and other parameters are consistent with those in the 3D model.
Generally speaking, the solute distribution around the static dendrite is relatively uniform, meaning the growth velocity of the dendrite trunks should be the same and should exhibit symmetrical growth in the six directions. Figure 4a demonstrates the streamline distribution as the 3D dendrite is falling; in the figure, the flow direction of the solution around the dendrite is upward. As the dendrite moves downward, the melt will wash the latent heat and solute released by the dendrite arm from the downstream tip to the upstream tip, meaning the growth of the dendrite is asymmetric.  In addition, as the dendrite falls, more solute is discharged from the upper tip of the dendrite and the diffusion speed of the discharged solute is not as fast as the accumulation speed of the solute, meaning there will be a high-concentration zone of solute enrichment in the upper part. The secondary dendrites also hinder the diffusion of the solute, resulting in the inhibition of the growth of the upper dendrite arm. This is consistent with the phenomenon observed by Badillo et al. [29].
As demonstrated in Figure 4a, in the 3D simulation process, the 3D streamline forms a circulation motion along the primary dendrite trunk in the vertical direction. In the horizontal direction, the streamline passes through the secondary dendrite arm, demonstrating the spatial fluid flow more realistically. In the 2D simulation process (as demonstrated in Figure 4b), a clockwise circulation motion is formed on the left side of the dendrite and a counterclockwise circulation motion symmetrical to the left side is formed on the right side of the dendrite. It can be observed that the solution is hindered when it flows to the dendrite trunk and can only bypass the dendrite tip along the dendrite arm, which makes In addition, as the dendrite falls, more solute is discharged from the upper tip of the dendrite and the diffusion speed of the discharged solute is not as fast as the accumulation speed of the solute, meaning there will be a high-concentration zone of solute enrichment in the upper part. The secondary dendrites also hinder the diffusion of the solute, resulting in the inhibition of the growth of the upper dendrite arm. This is consistent with the phenomenon observed by Badillo et al. [29].
As demonstrated in Figure 4a, in the 3D simulation process, the 3D streamline forms a circulation motion along the primary dendrite trunk in the vertical direction. In the horizontal direction, the streamline passes through the secondary dendrite arm, demonstrating the spatial fluid flow more realistically. In the 2D simulation process (as demonstrated in Figure 4b), a clockwise circulation motion is formed on the left side of the dendrite and a counterclockwise circulation motion symmetrical to the left side is formed on the right side of the dendrite. It can be observed that the solution is hindered when it flows to the dendrite trunk and can only bypass the dendrite tip along the dendrite arm, which makes it difficult for the solute to drain away, creating a high-concentration zone. Through this comparison, it is obvious that the calculation of dendrite movement and growth and solute transport under 3D conditions is more realistic. Figure 5 demonstrates the curves for the total solid fraction over time under 3D and 2D conditions, which represent the growth velocity of the dendrite. It can be observed that the growth velocity of the 3D equiaxed crystal is faster, which is about 2.7 times that of the 2D equiaxed crystal under the same conditions; this is mainly because the 3D model has more diffusion directions (along the Z axis) than the 2D model. Therefore, the solute discharged from the tip is more easily diffused, which makes the concentration gradient of the solute at the tip front larger than that of the 2D model. This is consistent with the conclusions drawn by Eshraghim et al. [30]. phenomenon observed by Badillo et al. [29].
As demonstrated in Figure 4a, in the 3D simulation process, the 3D streamline forms a circulation motion along the primary dendrite trunk in the vertical direction. In the horizontal direction, the streamline passes through the secondary dendrite arm, demonstrating the spatial fluid flow more realistically. In the 2D simulation process (as demonstrated in Figure 4b), a clockwise circulation motion is formed on the left side of the dendrite and a counterclockwise circulation motion symmetrical to the left side is formed on the right side of the dendrite. It can be observed that the solution is hindered when it flows to the dendrite trunk and can only bypass the dendrite tip along the dendrite arm, which makes it difficult for the solute to drain away, creating a high-concentration zone. Through this comparison, it is obvious that the calculation of dendrite movement and growth and solute transport under 3D conditions is more realistic. Figure 5 demonstrates the curves for the total solid fraction over time under 3D and 2D conditions, which represent the growth velocity of the dendrite. It can be observed that the growth velocity of the 3D equiaxed crystal is faster, which is about 2.7 times that of the 2D equiaxed crystal under the same conditions; this is mainly because the 3D model has more diffusion directions (along the Z axis) than the 2D model. Therefore, the solute discharged from the tip is more easily diffused, which makes the concentration gradient of the solute at the tip front larger than that of the 2D model. This is consistent with the conclusions drawn by Eshraghim et al. [30].  Figure 6 demonstrates a comparison of the falling velocity rates for 2D and 3D dendrites. It can be observed from the figure that the 3D dendrite and the 2D dendrite start to fall at the same time, but the falling velocity for the 3D dendrite is obviously higher than that of the 2D dendrite at the same time point. The main reason for this phenomenon is that the solute discharged from the 3D dendrite tip diffuses faster and more Cu is discharged from the dendrite tip. The density of Cu is higher than for Al, resulting in a downward flow field around the dendrite tip. The downward flow field around the 3D dendrite is stronger, causing the dendrite to fall faster. Figure 7 demonstrates the concentration distribution along the dendrite axis in 2D and 3D. In the figure, the ordinate represents the solute content, and the abscissa represents the grid number of the dendrite cross-section. The left side of the line represents the concentration distribution in front of the upstream tip and the right side represents the concentration distribution in the downstream tip. It can be observed from the figure that the 2D and 3D concentration distributions display a certain consistency. The lowest horizontal line in the figure represents the solute content in the solid phase inside the equiaxed crystal, with a low concentration of about 0.78%. The solute concentration increases sharply at the solid-liquid interface and the highest solute content reaches 5.9%. Then, the concentration in the liquid phase decreases with the distance from the dendrite. On the whole, the concentration at the front of the upstream tip is significantly lower than that at the downstream tip and the concentration gradient at the front of the upstream tip is greater than that at the downstream tip. It can be observed that the solute concentration at the front of the 3D dendritic tip is lower than that at the 2D dendritic tip, and the solute diffusion layer at the 3D dendritic tip is thinner than that at the 2D dendritic tip. This is because, in the process of 3D dendrite movement, the melt drives the solute at the upstream tip to flow to the downstream tip, and also accelerates the temperature diffusion, shortens the solute diffusion layer thickness at the interface front of upstream tip, and increases that at the front of downstream tip. Figure 6 demonstrates a comparison of the falling velocity rates for 2D and 3D dendrites. It can be observed from the figure that the 3D dendrite and the 2D dendrite start to fall at the same time, but the falling velocity for the 3D dendrite is obviously higher than that of the 2D dendrite at the same time point. The main reason for this phenomenon is that the solute discharged from the 3D dendrite tip diffuses faster and more Cu is discharged from the dendrite tip. The density of Cu is higher than for Al, resulting in a downward flow field around the dendrite tip. The downward flow field around the 3D dendrite is stronger, causing the dendrite to fall faster.

Figure 6.
Falling velocity rates for 2D and 3D dendrites. Figure 7 demonstrates the concentration distribution along the dendrite axis in 2D and 3D. In the figure, the ordinate represents the solute content, and the abscissa represents the grid number of the dendrite cross-section. The left side of the line represents the concentration distribution in front of the upstream tip and the right side represents the concentration distribution in the downstream tip. It can be observed from the figure that the 2D and 3D concentration distributions display a certain consistency. The lowest horizontal line in the figure represents the solute content in the solid phase inside the equiaxed crystal, with a low concentration of about 0.78%. The solute concentration increases sharply at the solid-liquid interface and the highest solute content reaches 5.9%. Then, the concentration in the liquid phase decreases with the distance from the dendrite. On the whole, the concentration at the front of the upstream tip is significantly lower than that at the downstream tip and the concentration gradient at the front of the upstream tip is greater than that at the downstream tip. It can be observed that the solute concentration at the front of the 3D dendritic tip is lower than that at the 2D dendritic tip, and the solute diffusion layer at the 3D dendritic tip is thinner than that at the 2D dendritic tip. This is because, in the process of 3D dendrite movement, the melt drives the solute at the upstream tip to flow to the downstream tip, and also accelerates the temperature diffusion, shortens the solute diffusion layer thickness at the interface front of upstream tip, and increases that at the front of downstream tip. The 3D dendrite arm is longer, which can affect the flow field in a larger rang it moves, meaning the solute discharged from the tip has a larger diffusion range. fluid can flow around the dendrite arm and wash the solute away directly from t of the dendrite arm, which greatly shortens the diffusion distance and reduces t centration enrichment.
In summary, there are many differences between the 2D model and the 3D m the calculation of the motion of equiaxed crystals. Obviously, the 3D CA-LBM mo more accurately describe the movement and growth of equiaxed crystals.

Study of the Movement Behavior of Multiple Dendrites
The calculation parameters in this section are consistent with those in Sect Figure 8 demonstrates the dendrite morphology and partial streamline distributio five dendrites in the stationary state and the moving state. It can be observed fr streamline distribution in Figure 8 that as the dendrite falls, due to the flow of the cluster of closed streamlines around the dendrite is formed. The streamline direc the 3D vortex are upward near the equiaxed crystal and downward far away fr equiaxed crystal, which means that the flow field always flows around the dend the moving state, the streamlines are distributed among the dendrites, which in that the interactions among the dendrites cause the solution convection. The indep fluid vortices begin to contact and gradually merge, the flow field in the same dire superimposed and enhanced, and the flow field in the opposite direction is partiall When multiple equiaxed crystals are falling, because there is no external forced The 3D dendrite arm is longer, which can affect the flow field in a larger range when it moves, meaning the solute discharged from the tip has a larger diffusion range. The 3D fluid can flow around the dendrite arm and wash the solute away directly from the side of the dendrite arm, which greatly shortens the diffusion distance and reduces the concentration enrichment.
In summary, there are many differences between the 2D model and the 3D model in the calculation of the motion of equiaxed crystals. Obviously, the 3D CA-LBM model can more accurately describe the movement and growth of equiaxed crystals.

Study of the Movement Behavior of Multiple Dendrites
The calculation parameters in this section are consistent with those in Section 4.1. Figure 8 demonstrates the dendrite morphology and partial streamline distribution of the five dendrites in the stationary state and the moving state. It can be observed from the streamline distribution in Figure 8 that as the dendrite falls, due to the flow of the melt, a cluster of closed streamlines around the dendrite is formed. The streamline directions of the 3D vortex are upward near the equiaxed crystal and downward far away from the equiaxed crystal, which means that the flow field always flows around the dendrite. In the moving state, the streamlines are distributed among the dendrites, which indicates that the interactions among the dendrites cause the solution convection. The independent fluid vortices begin to contact and gradually merge, the flow field in the same direction is superimposed and enhanced, and the flow field in the opposite direction is partially offset. When multiple equiaxed crystals are falling, because there is no external forced convection in this paper, the flow field produced by the natural convection driving the movement of dendrites is weak and multiple dendrites hinder each other, so the falling distance is small. that the interactions among the dendrites cause the solution convection. The independent fluid vortices begin to contact and gradually merge, the flow field in the same direction is superimposed and enhanced, and the flow field in the opposite direction is partially offset. When multiple equiaxed crystals are falling, because there is no external forced convection in this paper, the flow field produced by the natural convection driving the movement of dendrites is weak and multiple dendrites hinder each other, so the falling distance is small.   Figure 9 demonstrates the relationship between the movement velocity and time, which can quantitatively describe the movement process for the dendrites. The results indicate that all of the dendrites fall at first, then the acceleration decreases, and then the velocity decreases, whereby their falling velocity is less than that of a single dendrite. For dendrite 2, because it is far away from other dendrites, the interaction between dendrite 2 and other dendrites will be smaller, and the solute buoyancy and thermal buoyancy caused by the solute and heat discharged from its tip will have a greater impact on the falling velocity. For dendrite 5, it is close to dendrite 3. When dendrite 3 falls, the solute discharged from its tip will hinder the falling of dendrite 5 and the thermal buoyancy will cause an upward flow field around dendrite 5, which will slow down the falling velocity. For dendrite 3, because it is in the middle of these dendrites, it will be affected by the upward combined flow field of dendrites 1 and 2 and the downward flow field of dendrites 4 and 5; therefore, the dendrites will be affected by the 3D vortex generated by the movement of other dendrites in the process of falling, so the falling velocity curve presents a fluctuating trend.
Crystals 2021, 11, x FOR PEER REVIEW Figure 9 demonstrates the relationship between the movement velocity an which can quantitatively describe the movement process for the dendrites. The indicate that all of the dendrites fall at first, then the acceleration decreases, and t velocity decreases, whereby their falling velocity is less than that of a single dend dendrite 2, because it is far away from other dendrites, the interaction between de and other dendrites will be smaller, and the solute buoyancy and thermal bu caused by the solute and heat discharged from its tip will have a greater impac falling velocity. For dendrite 5, it is close to dendrite 3. When dendrite 3 falls, th discharged from its tip will hinder the falling of dendrite 5 and the thermal buoya cause an upward flow field around dendrite 5, which will slow down the falling v For dendrite 3, because it is in the middle of these dendrites, it will be affected upward combined flow field of dendrites 1 and 2 and the downward flow field drites 4 and 5; therefore, the dendrites will be affected by the 3D vortex generated movement of other dendrites in the process of falling, so the falling velocity curve p a fluctuating trend. The 3D vortex around a dendrite hinders its vertical drop, and because the str around the dendrite not only flows vertically down along the Z axis but also flows the dendrite, the surrounding solution will produce lateral force on the dendrit horizontal plane and the dendrite will also move in the horizontal direction.
The lateral displacement of the dendrites is presented in Tables 2 and 3.  The 3D vortex around a dendrite hinders its vertical drop, and because the streamline around the dendrite not only flows vertically down along the Z axis but also flows around the dendrite, the surrounding solution will produce lateral force on the dendrite in the horizontal plane and the dendrite will also move in the horizontal direction.
The lateral displacement of the dendrites is presented in Tables 2 and 3.
Here, − means that the dendrite moves to the negative direction of the coordinate axis and + means that the dendrite moves to the positive direction of the coordinate axis. The displacement distance is calculated by the position of the center of mass. According to the displacement of dendrites, it can be observed that the lateral forces acting on the dendrites at different positions are not the same. Under the action of the flow field generated by other dendrites, dendrite 1 will move to the positive directions of the x-axis and y-axis, while dendrite 2 will move to the positive direction of the x-axis and negative direction of the y-axis. Compared with dendrite 1, dendrite 2 is in a more central position and is closer to other dendrites, so it has a stronger effect and longer horizontal displacement distance. Other dendrites are similarly affected, so will not be discussed here.

Conclusions
In this paper, the LBM-CA model, Ladd method, non-equilibrium extrapolation method, and Newton's law of classical mechanics were combined to establish a 3D model for calculating the motion and growth process of equiaxed crystals for the first time. The model was verified and the motion of single and multiple 3D equiaxed crystals was calculated. Compared with the 2D method, the difference between the two growth models was analyzed. The results demonstrated that: (1) In the 3D simulation, the 3D eddy current drove the solute to diffuse around the dendrite, so the solute diffusion was faster, the solute concentration between dendrites was lower, and the falling and growth speeds were faster than in the 2D model. Therefore, the 3D model more truly reflected the movement of equiaxed crystals and the melt flow around dendrites. (2) The growth modes of dendrites in motion are different from that in static state, and multiple dendrites will influence each other in the processes of movement and growth, changing the original movement and growth modes.

Conflicts of Interest:
The authors declare no conflict of interest.