Numerical Simulation of Dendrite Motion Fidelity Based on the Interface Capture Method

,


Introduction
Te solidifcation of a large casting product has an important impact on its fnal quality, and the microstructure obtained by solidifcation directly afects the performance of the material [1]. During the solidifcation of a large casting product, many grains undergo various behaviors, such as moving, colliding, separating, sinking, and accumulation, inside the melt [2,3]. Many free grains (mainly equiaxed grains) constantly move in the melt under natural convection or forced convection, and the motion behaviors of these free grains have important efects on the fnal microstructure, tissue distribution, and component distribution of the product, thereby ultimately afecting its performance [4][5][6].
To date, the numerical simulation methods of microstructures and dendrite growth processes in solidifcation mainly include phase feld (PF) and cellular automata (CA) methods [7]. Most of the simulations for dendrite motion use the PF method because the dendrite interface calculated by the PF is a continuous interface, and it is easier to trace the solid-liquid boundary when the dendrite is moving. Rojas et al. [8] successfully simulated the motion process of a single dendrite by coupling PF with a lattice Boltzmann method (LBM) and equations of motion. Qi et al. [9] calculated the growth and motion processes of multiple dendrites by combining the extended phase feld model with the N-S equation; however, this model does not address the collision behaviors between multiple dendrites. Sakane et al. [10][11][12][13] used multi-GPU parallel computing to greatly improve the scale and efciency of computing, and they simulated the solidifcation processes of hundreds of equiaxed crystal motion collision behaviors. However, since the PF method relies on a small mesh size, the calculation amount is very large and the calculation efciency is low when calculating the macroscopic segregation of a large-scale ingot.
Te CA method is simple and fast to calculate, and it is easy to combine with various physical processes, which has broad industrial prospects; however, there is less research on dendrite motion, and it is still in its infancy. Liu et al. [14] and Wu et al. [15] calculated the falling behavior of a single dendrite during solidifcation using the CA-LB method coupled with the Ladd method. Te solid-liquid interface calculated by the CA model is sharp; when the dendrite motion is calculated using the CA model of Liu et al., the movement of the solid-liquid interface in space moves on a grid scale, and it is not continuously calculated by time steps. When calculating the horizontal and vertical movement characteristics of the dendrite with this method, the shape of the dendrite is maintained; however, in the calculation of oblique movement, especially when calculating the circular motion, with increasing calculation time, the dendrite undergoes signifcant deformation, the morphology of the dendrite is not maintained, and the solute accumulates at the interface after a certain period of calculation, afecting the growth morphology of the dendrite. To realize the continuous movement of the interface according to the time step, Zhang et al. [16] developed dynamic grid technology to process the movement of dendrites. Tis method maintains the morphologies of dendrites well. However, when addressing multiple dendrites, the dynamic grid scheme greatly increases the calculation amount and reduces the calculation efciency.
In this paper, the numerical model level set method and volume of fuid (VOF) method of free interface tracking in computational fuid mechanics are used to track the solidliquid interface during the dendrite movement; the fdelity degree of the dendrite movement morphology calculation using these two methods is investigated. Te two methods are coupled with the CA-LB model to simulate the movement of dendrites under convection. Te coupled model is used to track the solid-liquid interface of moving dendrites with neglected growth, and the optimal interface tracking method is selected by analyzing the morphology loss rate.

Model Description
In this paper, the CA method is used to calculate the growth of equiaxed crystals and the LB model is used to calculate the transport of the temperature feld, fow feld, and solute feld in the melt where the equiaxed crystals are located. Te specifc algorithm is detailed in the literature [16]. Ten, the displacements of dendrites are obtained by solving the equation of motion, and the position of the solid-liquid interface is updated by using the level set method, SLIC-VOF method, and PLIC-VOF method.

Level Set Method. Te level set method was proposed by
Osher and Sethian in 1988 to solve the motion of two-phase fow [17]. Te basic idea is to defne the level set function ϕ(x, t) as the distance function from the current node to the interface, and the interface is the zero level of the distance function. Te zero level at the new moment is calculated by solving the development equation of the level set function to track the interface development and change.
Te distance function ϕ(x, t) is defned as follows: where d(x, Γ) is the shortest distance between the point and the interface, Γ represents the solid-liquid interface, and Ω 1 and Ω 2 represent the liquid and solid phases, respectively. In the study of two-phase fow, the change in the level set function follows the convection transport equation: Te motion of the interface is tracked by solving equation 2. Te position of the interface is realized by fnding the set of points that always satisfy φ � 0. Te direction and curvature of the outer normal of the free interface are expressed as follows: However, in the process of solving the level set equation, due to the inherent dissipation efect of the numerical method, when the contour line of φ � 0 is calculated with several time steps, it no longer meets the defnition of equation (1), resulting in the deformation of the solved moving interface [18]. Terefore, the level set equation needs to be constantly reinitialized.
When set at time t, the level set function φ 0 is obtained. To reconstruct the function ϕ(x, t), two conditions must be satisfed: (1) ϕ(x, t) satisfes equation (1), which is the sign distance function (2) ϕ(x, t) and φ 0 have the same zero level Tese two conditions are satisfed by solving the stable solution of the following initial value problem: where sign(φ 0 ) is a symbolic function, defned as follows [19]: where ε is a minimum value and the denominator cannot equal zero. Te distance function under the new zero level is obtained by iterating equation (4) several times until it is stable.

VOF Method.
Te VOF model is an interface tracking method proposed by Hirt and Nichols [20] in 1981. Tis model uses the concept of the fuid volume function F to represent the proportions of diferent phases in a unit, and it restructures the interface through the change in the fuid volume function over time to determine the position, shape, and change in the moving interface at every moment for interface tracking. In this article, F is used to represent the volume fraction of the solid phase in the grid. Te transfer equation of the volume function F is calculated by the following formula: Te key to tracking the interface using the VOF method is to reconstruct the free interface; that is, the specifc position of the interface in the mesh according to the volume function F of the mesh cell and its neighboring mesh cells is reconstructed. To date, there are many methods for VOF interface reconstruction, such as the SLIC method [20], FLAIR method [21], PLIC method [22], and CICSAM method [23]. In this paper, the SLIC method and PLIC method are adopted to track the motion interface, and the accuracies of the two methods are compared. Te following introduces the two interface reconstruction methods.

SLIC Method.
Te SLIC method defnes the free interface in the grid as horizontal and vertical and treats the fuid-free interface as a local single value function Y(x) and X(y). Te nine-grid template (as shown in Figure 1) is used to calculate the Y l value of grid columns i − 1, i + 1, and the X l value of grid columns j − 1, j + 1 through equations (7) and (8). Te slope values dY/dx and dX/dy of the free surface on each grid are calculated by equation (9), and then the position and direction of the free interface on grid (i, j) are determined according to the volume function and the magnitude of the slope.
By comparing the absolute value of dY/dx and dX/dy size, the interface in the grid is determined to be horizontal or vertical. If | (dY/dx)| < | (dX/dy)|, the interface is considered horizontal; otherwise, the interface is vertical.

PLIC Method.
Te basic principle of the PLIC method is as follows. In a single grid, a frst-order line is used to approximate the phase interface; that is, the dip angle between the moving interface and the boundary line is determined. Te slope and position of the line are determined by using the dip angle and the volume function in the grid, and the phase interface in the grid is constructed. Ten, the fuid volume function values of the central grid and the neighboring grids are modifed by calculating the volume fux of the fuid fowing through the central grid boundary to the neighboring grids after a time step.
If the cell is a full grid (F � 1), it is frst determined whether it is the donor or the recipient by the velocity direction of the four boundaries. Ten, the volume of fuid fowing through adjacent grids around the boundary is calculated in a time step. When the cell is a half grid (0 < F < 1), the normal vector n � (n x i,j , n y i,j ) of the interface in the grid is calculated frst. Te calculation formula is as follows: Ten, according to the normal vector, the angle β between the motion interface and the X-axis is determined as follows: Te angle α is defned as follows: In this manner, 20 diferent combinations of the interface in the grid are obtained. Except for the four simple cases when n x i,j � 0 and n y i,j � 0, the other 16 kinds of symmetry and inversion are simplifed to four types (see Figure 2). Tis angle and the volume function are used to determine which type of symmetry or inversion is apparent; then, the slope and position of the line are calculated, the interface within the mesh is constructed, and the volume of fuid fowing through adjacent meshes within a time step is calculated across the surrounding boundaries. Finally, the fuid volume function of this grid and adjacent grids is modifed and solved iteratively.  Figure 3(a). When calculating the rotating feld, the initial position of the square center is at [0.5, 0.5], and the rotation velocity feld is (u, v) � (− π(y − 0.5), π(x − 0.5)), as shown in Figure 3

Results and Discussion
Te solution error is quantifed as follows: where F n i,j is the calculated solution after n time steps and F 0 i,j is the initial solution. Te calculation errors of the abovementioned three methods in the two fow felds are shown in Table 1. Figure 4 shows the simulation results of the translational feld. Te interface calculated by the level set method is relatively clear; however, the sharp corners of the interface appear smooth, the edges gradually disappear, and the deformation is serious. Figure 5 shows the mass change curve of the solid square motion process. Te fgure shows that with increasing numbers of calculation steps, the block mass calculated by the level set method decreases continuously, and the mass loss is serious.
Te sharpness of the interface calculated by the SLIC method is the best, and the interface is clear. Figure 5 shows that the method has a small loss rate of morphology and a small mass change range. Table 1 shows that the numerical dissipation calculated by the PLIC method is the smallest at only − 0.0004%. However, the PLIC method needs to calculate the normal of the interface at the moving interface unit when reconstructing the interface; thus, there is a slight smoothing at the sharp corners of the interface.
Te simulation results of the rotating feld are shown in Figure 6 for each of the three methods. As shown in Figures 6(a)-6(c), when the level set method is used to process the rotating feld, the smooth interface phenomenon is more serious, the interface edges and corners are completely smoothed, and the area loss is large. In addition, when using the level set method to calculate the interface movement, it is necessary to determine the specifc location of the solid-liquid boundary in the interface unit to initialize the distance function. However, the specifc location of the dendrite boundary in the unit calculated by the CA method is unknown. In the actual calculation, the growth and movement process of dendrites are carried out simultaneously, and the profles of dendrites change at each time step. If the level set method is used to calculate the dendrite movement, initialization is required at each time step, greatly increasing the calculation time. By considering the above reasons, the level set method is difcult to couple with the CA model and is not suitable for calculating the movements of dendrites. Figures 6(d)-6(f ) show that the SLIC method has a large error when calculating the rotating feld, the fneness of the graphic surface worsens, the corner is smoothed, and the volume expands obviously. Tis phenomenon occurs because the method has zero order precision and only handles simple motion processes.
As seen from the calculation results of Figures 6(g)-6(i), the PLIC method maintains the original morphology well. Although there is slight smoothing at the corner, the overall fneness is signifcantly higher than that of the previous two interface tracking methods. By comparing the morphology loss rate in Table 1, the error of the PLIC method is only − 0.000 4%, which is much smaller than that of the other two methods; this error does not have a signifcant impact on the subsequent calculation.

Simulation of the Dendrite Motion Process.
As seen from the calculation results in Section 3.1, the error of the SLIC method is small in the calculation of the translational feld, and the PLIC method obtains relatively high precision in the calculation of the translational feld and rotating feld. However, due to the complex interface morphologies of dendrites, to verify the tracking abilities of the above-given two VOF methods for complex moving interfaces, the above-given two methods are coupled with the CA-LB model to calculate the motion behaviors of dendrites with negligible growth in translational and rotational felds. Te whole calculation domain is [1, 0] × [1, 0], the grid is 400 × 400, and the calculation time step is δt � 0.0005s. To facilitate the research, the initial fow feld is set as a static state, the initial undercooling is 7 K, the boundary condition of the temperature feld is an adiabatic boundary condition, and the solute feld is a nondifusion boundary condition. In the calculation of the translational feld, the initial nucleation point is at [0.25, 0.75], and the preferred growth angle is 0°. After 1.5 s of dendrite growth, the dendrites stop growing and are given an initial velocity of (u, v) � (0.4, − 0.4). Te initial position and morphology before movement are shown in Figure 7(a).
Te SLIC method and PLIC method are used to trace the movement interfaces of dendrites. Figures 8(a)-8(c) are the calculation results of the SLIC method, and the dendrites are seriously deformed and lose their original morphological characteristics. Figure 9 shows that the dendrite mass is seriously expanded. Tis phenomenon occurs because the interface reconstructed by the SLIC method is rough. When the calculated interface shape is complex, the reconstructed interface is quite diferent from the actual interface; the actual volume function transmission cannot be accurately calculated, resulting in serious deformation. Figures 8(d)-8(f ), the PLIC method has good adaptability to complex interfaces. Te calculated dendrite morphology is maintained well, and only slight smoothing occurs at the tip of the dendrite. As shown in Table 2, the morphology loss rate is only − 0.001%.

As shown in
In the calculation of the rotating feld, the initial nucleation point is at [0.25, 0.75], and the preferred growth angle is 0°. After 1.5 s of dendrite growth, the dendrites stop growing and are given an initial velocity of (u, v) � (− π(y − 0.5), π(x − 0.5)). Te other conditions are the same as the translational feld calculations. Te initial morphologies of the dendrites are shown in Figure 7(b). Figures 10(a)-10(f ) show the morphologies and positions of dendrites after rotations of 2π, 4π, and 8π calculated by the SLIC and PLIC methods, respectively. Te fgure shows that the rotation of the complex interface calculated by the SLIC method is the same as that of the simple interface. After rotation, the interface is seriously blurred, the tip of the dendrite is smoothed, the volume is seriously expanded, and the direction of the dendrite is obviously shifted. Te calculation results of the PLIC method maintain the original morphology, and the dendrite tip appears slightly distorted in the initial stage of rotation; however, with the progress of rotation, this deformation tends to be stable. According to Table 2, the appearance loss rate is only − 0.001%, and the error does not cause additional interference to the calculation. Te PLIC method is similar to the CA method in that it calculates based on the solid fraction within the grid, and the solidliquid interface obtained by the PLIC method is sharp, so it does not afect the calculation of dendrite growth when coupled with the CA method. When calculating the movement of dendrites including growth, at the same time   step, after calculating the growth of dendrites using the CA method, the position of the dendrites after movement can be calculated using the PLIC method based on the velocity of dendrite movement and the solid fraction in the grid, realizing the coupling of the CA and PLIC methods. By coupling the PLIC method with the CA-LB model, the dendrite motion behavior can be continuously calculated under a single set of grids.

Dendrite Recombination Process.
Te movement process of grains during ingot forming is complex and changeable.
To verify the fdelity of the PLIC method in simulating the real movement of grains, the rotational and translational processes are combined to calculate the grain movement.
Te initial condition setting of the computing domain is the same as in Section 3.2. Te initial nucleation point of the dendrite is located at [0.5, 0.75]. After 1.5 s of dendrite growth, the dendrites stop growing and have an initial angular velocity of ω � 2π, giving the whole calculation domain (u, v) � (− π(y − 0.5), π(x − 0.5)) velocity feld. Te dendrite is allowed to rotate around the center of the computational domain while rotating. After 4 s of movement, the dendrites return to the initial position; the movement process is shown in Figure 11(a). As shown in Figure 11(a), the dendrite tip gradually becomes smooth, and the tip radius increases in the process of movement until it reaches a stable state. Tis phenomenon is because the interface normal of the cell must be solved in the interface reconstruction of the PLIC method; however, the interface curvature of the dendrite tip is large, and Formula (9) has only frst-order accuracy. Terefore, there are errors in the interface reconstruction of the dendrite tip cells, resulting in a smooth interface. Figure 11(b) shows the comparison of the dendrite contour before and after the movement; the blue solid line is the dendrite contour before the movement, and the red dotted line is the dendrite contour after the movement. After the end of the movement, the position and direction of the dendrites do not shift, the dendrite tip has slight smoothing, and the overall contour before and after the movement is not much diferent. Tis motion process combines most of the motion behaviors that may exist in the real solidifcation of dendrites. Te calculation results show that the PLIC method calculates the motion behaviors of dendrites with small errors and solve the motion discontinuity phenomenon in the previous motion model.

Conclusion
In this paper, three commonly used interface capture methods are selected to track the moving interface in the translational feld and the rotating feld, and the interface tracking ability and error size of diferent models are compared. Te calculation results are as follows: (1) Te level set method is not suitable for calculating the motion behaviors of dendrites because of its weak tracking ability of the moving interface of rigid objects and its difculty in coupling with the CA model; the SLIC method is rough and only calculates objects with simple interfaces, and the error in the calculation of the rotating feld is large. Te PLIC method is better than the previous two methods, and the interface is more refned. Although there is a slight smoothing at sharp corners, the overall error is small. In addition, the VOF method is more convenient for coupling with the CA model, which is suitable for simulating dendrite motion behavior.
(2) Among the VOF methods, the calculation results of the oblique translation and rotational movements of the equiaxed dendrites calculated by the SLIC and PLIC methods, respectively, show that the deformation of the equiaxed dendrites calculated by the SLIC method is relatively serious and the error is large; the equiaxed dendrites calculated by the PLIC method have good fdelity, regardless of oblique movement or rotation, and the morphology loss rate reaches − 0.001%.

Data Availability
Te data used to support the fndings of this study are available from the corresponding author upon request.

Conflicts of Interest
Te authors declare that they have no conficts of interest. Journal of Engineering 9