Abstract
During alloy solidification, many free grains in the melt have important effects on the final microstructure and composition distributions. In this paper, grain motion is calculated based on an interface tracking method coupled with a cellular automata (CA) method. First, the interface tracking capabilities of the level set, simple linear interface calculation (SLIC), and piecewise linear interface calculation (PLIC) methods are compared, and the fidelity of the three models is explored. Then, the coupling degrees of these three models with the CA method are analyzed. Finally, the PLIC method is used to simulate various behaviors of grain movement and to verify the authenticity of the dendrite motion calculation. The simulation results show that the VOF methods more readily coupled with the CA model than the level set method, and it is more suitable for calculating the motion behaviors of dendrites. Among the VOF methods, the interface reconstructed by the SLIC method is relatively rough and can only calculate objects with simple morphologies. The PLIC method has a fine interface and small error in the calculation of dendrite movement, and it does not significantly impact the subsequent calculations.
1. Introduction
The solidification of a large casting product has an important impact on its final quality, and the microstructure obtained by solidification directly affects the performance of the material [1]. During the solidification 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 effects on the final microstructure, tissue distribution, and component distribution of the product, thereby ultimately affecting its performance [4–6].
To date, the numerical simulation methods of microstructures and dendrite growth processes in solidification mainly include phase field (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 field model with the N-S equation; however, this model does not address the collision behaviors between multiple dendrites. Sakane et al. [10–13] used multi-GPU parallel computing to greatly improve the scale and efficiency of computing, and they simulated the solidification 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 efficiency is low when calculating the macroscopic segregation of a large-scale ingot.
The 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 solidification using the CA-LB method coupled with the Ladd method. The 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 significant deformation, the morphology of the dendrite is not maintained, and the solute accumulates at the interface after a certain period of calculation, affecting 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. This 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 efficiency.
In this paper, the numerical model level set method and volume of fluid (VOF) method of free interface tracking in computational fluid mechanics are used to track the solid-liquid interface during the dendrite movement; the fidelity degree of the dendrite movement morphology calculation using these two methods is investigated. The two methods are coupled with the CA-LB model to simulate the movement of dendrites under convection. The 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.
2. 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 field, flow field, and solute field in the melt where the equiaxed crystals are located. The specific algorithm is detailed in the literature [16]. Then, 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.
2.1. Level Set Method
The level set method was proposed by Osher and Sethian in 1988 to solve the motion of two-phase flow [17]. The basic idea is to define the level set function as the distance function from the current node to the interface, and the interface is the zero level of the distance function. The 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.
The distance function is defined as follows:where is the shortest distance between the point and the interface, represents the solid-liquid interface, and and represent the liquid and solid phases, respectively.
In the study of two-phase flow, the change in the level set function follows the convection transport equation:
The motion of the interface is tracked by solving equation 2. The position of the interface is realized by finding the set of points that always satisfy . The 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 effect of the numerical method, when the contour line of is calculated with several time steps, it no longer meets the definition of equation (1), resulting in the deformation of the solved moving interface [18]. Therefore, the level set equation needs to be constantly reinitialized.
When set at time t, the level set function is obtained. To reconstruct the function , two conditions must be satisfied:(1) satisfies equation (1), which is the sign distance function(2) and have the same zero level
These two conditions are satisfied by solving the stable solution of the following initial value problem:where is a symbolic function, defined as follows [19]:where ε is a minimum value and the denominator cannot equal zero. The distance function under the new zero level is obtained by iterating equation (4) several times until it is stable.
2.2. VOF Method
The VOF model is an interface tracking method proposed by Hirt and Nichols [20] in 1981. This model uses the concept of the fluid volume function F to represent the proportions of different phases in a unit, and it restructures the interface through the change in the fluid 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. The transfer equation of the volume function F is calculated by the following formula:
The key to tracking the interface using the VOF method is to reconstruct the free interface; that is, the specific 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. The following introduces the two interface reconstruction methods.
2.2.1. SLIC Method
The SLIC method defines the free interface in the grid as horizontal and vertical and treats the fluid-free interface as a local single value function and . The nine-grid template (as shown in Figure 1) is used to calculate the value of grid columns i − 1, i + 1, and the value of grid columns j − 1, j + 1 through equations (7) and (8). The slope values and 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 and size, the interface in the grid is determined to be horizontal or vertical. If , the interface is considered horizontal; otherwise, the interface is vertical.
2.2.2. PLIC Method
The basic principle of the PLIC method is as follows. In a single grid, a first-order line is used to approximate the phase interface; that is, the dip angle between the moving interface and the boundary line is determined. The 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. Then, the fluid volume function values of the central grid and the neighboring grids are modified by calculating the volume flux of the fluid flowing through the central grid boundary to the neighboring grids after a time step.
If the cell is a full grid (F = 1), it is first determined whether it is the donor or the recipient by the velocity direction of the four boundaries. Then, the volume of fluid flowing 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 of the interface in the grid is calculated first. The calculation formula is as follows:
Then, according to the normal vector, the angle β between the motion interface and the X-axis is determined as follows:
The angle is defined as follows:
In this manner, 20 different combinations of the interface in the grid are obtained. Except for the four simple cases when and , the other 16 kinds of symmetry and inversion are simplified to four types (see Figure 2). This 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 fluid flowing through adjacent meshes within a time step is calculated across the surrounding boundaries. Finally, the fluid volume function of this grid and adjacent grids is modified and solved iteratively.
3. Results and Discussion
3.1. Simulation of the Solid Square Motion Process
To compare the abilities of the abovementioned models to track the moving interface, the level set method, SLIC method, and PLIC method are used to analyze and compare the motion deformation of the solid square in the constant flow field and the rotating flow field. By comparing with the initial morphology and calculating the mass loss, the tracking capabilities of the abovementioned three models are evaluated.
The whole calculation domain is [1, 0] × [1, 0], the grid is 400 × 400, the calculation time step is , and the initial square size is 0.125 × 0.125. For the level set method, the boundary line of the square is a zero contour line, the inside of the square is the region, and the outside is the region. For the VOF class method, the square internal volume function F is assigned a value of 1.0 and an outer value of 0.0. During the calculation of the translational flow field, the initial position of the square center is at [0.25, 0.75], and the constant velocity field is , as shown in Figure 3(a). When calculating the rotating field, the initial position of the square center is at [0.5, 0.5], and the rotation velocity field is , as shown in Figure 3(b).
The solution error is quantified as follows:where is the calculated solution after n time steps and is the initial solution. The calculation errors of the abovementioned three methods in the two flow fields are shown in Table 1.
Figure 4 shows the simulation results of the translational field. The 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. The figure 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.
The 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.
The simulation results of the rotating field 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 field, 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 specific location of the solid-liquid boundary in the interface unit to initialize the distance function. However, the specific 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 profiles 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 difficult 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 field, the fineness of the graphic surface worsens, the corner is smoothed, and the volume expands obviously. This 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 fineness is significantly 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 significant impact on the subsequent calculation.
3.2. 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 field, and the PLIC method obtains relatively high precision in the calculation of the translational field and rotating field. 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 fields.
The whole calculation domain is [1, 0] × [1, 0], the grid is 400 × 400, and the calculation time step is . To facilitate the research, the initial flow field is set as a static state, the initial undercooling is 7 K, the boundary condition of the temperature field is an adiabatic boundary condition, and the solute field is a nondiffusion boundary condition. In the calculation of the translational field, 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 . The initial position and morphology before movement are shown in Figure 7(a).
The 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. This phenomenon occurs because the interface reconstructed by the SLIC method is rough. When the calculated interface shape is complex, the reconstructed interface is quite different from the actual interface; the actual volume function transmission cannot be accurately calculated, resulting in serious deformation.
As shown in Figures 8(d)–8(f), the PLIC method has good adaptability to complex interfaces. The 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%.
In the calculation of the rotating field, 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 . The other conditions are the same as the translational field calculations. The 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. The figure 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. The 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. The PLIC method is similar to the CA method in that it calculates based on the solid fraction within the grid, and the solid-liquid interface obtained by the PLIC method is sharp, so it does not affect 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.
3.3. Dendrite Recombination Process
The movement process of grains during ingot forming is complex and changeable. To verify the fidelity of the PLIC method in simulating the real movement of grains, the rotational and translational processes are combined to calculate the grain movement. The initial condition setting of the computing domain is the same as in Section 3.2. The 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 velocity field. The 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. This 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 first-order accuracy. Therefore, 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 different. This motion process combines most of the motion behaviors that may exist in the real solidification of dendrites. The 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.
4. Conclusion
In this paper, three commonly used interface capture methods are selected to track the moving interface in the translational field and the rotating field, and the interface tracking ability and error size of different models are compared. The calculation results are as follows:(1)The 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 difficulty 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 field is large. The PLIC method is better than the previous two methods, and the interface is more refined. 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 fidelity, regardless of oblique movement or rotation, and the morphology loss rate reaches −0.001%.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was funded by National Natural Science Foundation of China (Grant no. 51975182).