Abstract

In order to simulate the ice crystal growth with five-fold symmetry, a new phase-field model was established. The finite difference method was used to solve the governing equation, reproduce the growth morphology of ice crystals with five-fold symmetry, and study the influence of interfacial energy anisotropy and undercooling on the growth mechanism of ice crystals. The results indicate that when the strength of interfacial energy anisotropy () is greater than the critical value (1/24), the interface discontinuity and corners form at the tips of the main stem and side branches. The characteristic parameters of the ice crystal tip in the horizontal direction have an extreme value at equal to 0.1. With an increase in undercooling, the main branch becomes slender and side branches appear, which is consistent with the mechanism of crystal continuous growth. When undercooling reaches 0.4, ice crystal growth changes from diffusion control to dynamic control. Therefore, by selecting the quick freezing conditions, undercooling is controlled at around 0.4 and the interfacial energy anisotropy at around 0.1, and the tip radius of curvature is minimized so as to prevent large ice crystals from damaging the cell wall and leading to the decline in the quality of quick-frozen food. The findings of this study can help us have a better understanding of ice crystal growth in the process of quick-frozen food preservation.

1. Introduction

In the process of quick-frozen food preservation, the liquid water in cells condenses into ice crystals. Large ice crystals [13] destruct the cellular structure, which leads to a decline in food quality. Different shapes of ice crystals have different degrees of destruction for the cellular structure. Therefore, a study on the growth behavior of ice crystals in the process of quick-frozen food preservation has important practical significance for improving the quality of quick-frozen food.

The phase-field method and first principles are important methods for microstructure simulation of solidification. Zin et al. [4] used the first principles to determine the structure of magnesium borohydride. The phase-field method successfully avoids tracking the complex solid-liquid interface by introducing order parameters, and the energy coefficient is used to describe the anisotropy of the solid-liquid interface, which has obvious advantages in the study of the solid-liquid phase transition. Under the condition of weak interfacial energy anisotropy, all orientations appear on the equilibrium shape and crystal morphology is smooth. With an increase in interfacial energy anisotropy, some crystal orientations disappear from the equilibrium shape. In the study of metal solidification, the phase-field method was used to simulate the solidification microstructure formation process of Cu70Ni30 and Fe-C alloys [57], respectively. Ferreira et al. [8] obtained the two-dimensional morphology of Al-Cu alloys by the phase-field method, which is consistent with the experimental results. Qi et al. [9] used the phase-field crystal method to simulate the low-angle asymmetric tilt grain boundary structure and dislocation motion on a nanoscale. Yuan and Ding [10, 11] studied the faceted dendrite growth behavior of Ni-Cu alloys under forced flow. The effects of interfacial energy anisotropy strength and flow velocity on the growth of faceted dendrites in pure materials and alloys were evaluated [1214]. Yuan et al. [15, 16] simulated the dendrite growth behavior of HCP materials with strong interfacial energy anisotropy of the phase-field model.

In terms of ice crystal growth, Li et al. [17] reproduced the ice crystal growth process in supercooled solutions by the phase-field model for the first time. Chen et al. [18, 19] applied the phase-field model to the ice crystal growth of freeze concentration. Based on the law of ice crystal growth obtained from the simulation, we initially determined the method to reduce the solute entrainment in ice crystals. Han et al. [20, 21] applied the six-fold symmetrical phase-field model to explore the growth process of ice crystals on the surface of polar ships and revealed the formation mechanism of ice crystals on the surface of supercooled solids. Liu et al. [22] proposed that ice crystals were mainly in four-fold and six-fold symmetrical structures and that other symmetrical structures also existed. Yang et al. [23, 24] constructed a five-fold symmetrical phase-field model based on the Wheeler model to explore the effect of freezing time on ice crystal morphology, and the simulation results were consistent with the experimental results.

In this paper, based on the Eggleston method and the Karma model, five-fold symmetry of ice crystal growth modified to the interfacial energy anisotropy function was established, and a detailed correction process was given. The modified phase-field model is used to primarily study the growth behavior of ice crystals when the strength of interfacial energy anisotropy is greater than the critical value. The values of the strength of interfacial energy anisotropy and the undercooling degree with the smallest radius of curvature of ice crystals are summarized, which can provide reference for improving the technology of quick freezing and preservation of food.

The general sketch of the ice crystals with five-fold symmetry by the phase-field method is shown in Figure 1. First, the model is established. Then, the model is solved by the finite difference method and simulated on the computer server with the c++ language programming. Next, data are visualized by using Tecplot and Origin software. Finally, the data are analyzed to summarize ice crystal growth rules, and the optimal process parameters are obtained.

2. Interfacial Energy Anisotropy and Equilibrium Shape

In ice crystal growth, the energy gradient coefficient is used to reflect the effect of interfacial energy anisotropy, which is expressed as follows:where θ is the angle between the X-axis and the normal direction of the dendrite solid-liquid interface, is the constant, f(θ) is the energy gradient coefficient, and is the strength of interfacial energy anisotropy for five-fold symmetry.

2.1. Function of Interfacial Energy Anisotropy

Figure 2 shows the relationship between f(θ) and θ with different in rectangular coordinates. It can be concluded from plots that the curve of f(θ) produces an ever-repeating wave of peaks and that valleys have a period of 2/5π. When θ is 0°, 72°, 144°, 216°, and 288°, it is the peak that corresponds to the maximum value, and when θ is 36°, 108°, 180°, 252°, and 324°, it is the valley that corresponds to the minimum value. The peaks and valleys are 1.04 and 0.96, respectively, when is 0.04, and the peaks and valleys are 1.1 and 0.9, respectively, when is 0.1.

Figure 3 shows missing orientations using the 1/f(θ) plot and tangent lines when is 0.1. According to the inversion method [25] related to curve correction, if the 1/f(θ) plot is nonconvex, then its tangent lines will lie inside, and orientations will disappear from the equilibrium shape. If the 1/f(θ) plot is convex, then its tangent lines will lie outside, and orientations will appear on the equilibrium shape. For the nonconvex part of the 1/f(θ) plot, its convexity can be changed by adding tangent lines. The portions of the 1/f(θ) plot between these tangent lines and the origin will correspond to missing orientations, and all other orientations will appear in crystal morphology.

When the 1/f(θ) plot is nonconvex, interface stiffness functions S(θ) satisfy

Figure 4 shows the shape curves of S(θ) under different ε5 in polar coordinates. When ε5 is 0.04, all the values of S(θ) are positive, and the curve shape displays typical five-fold symmetry (see Figure 4(a)). When ε5 > 1/24, the values of S(θ) are negative at some positions, and the nonconvex part of the curve forms a convex part on the other side, corresponding to missing orientations (see Figure 4(b)).

Restricting attention to , missing orientations occur in the range of , and all orientations grow normally in the range of . In order to obtain accurate equilibrium crystal morphology, we adopt the method reported by Eggleston to regularize the function of interfacial energy anisotropy as follows [26]:

Regularized interface stiffness is as follows:

Critical missing orientations satisfy

We can get the value of by solving Equation (5).

2.2. Equilibrium Shape

The interfacial energy anisotropy function is substituted into the parameter equation of the equilibrium shape of crystals [27]; we obtain the following formulas:where A is a constant and X and Y are dimensionless distances.

The equilibrium shape of crystals under different is shown in Figure 5. It can be concluded from plots that, when ε5 ≤ 1/24, the equilibrium shape is smooth and continuous, and the overall shape is quasi-pentagonal. When ε5 > 1/24, some crystal orientations near the vertex of the pentagon disappear from the equilibrium shape, and the interface becomes discontinuous, which causes distortion “ears.” With an increase in , the range of crystal orientation disappearance increases, which resulted in an enlargement of these “ear” parts. After is regularized, “ear” parts degenerate into a point.

3. Phase-Field Model

The Eggleston method was used to modify the interfacial energy anisotropy function in the Karma model; a new phase-field model was established for ice crystal growth with five-fold symmetry. The expression of the governing equation is given for as follows:For , we obtainwhere the phase-field parameters are the same as in [23].

4. Numerical Simulations

An initial crystal is in the center of size 1200 × 1200 grid computational domains, and its radius is assumed to be R0. The initial conditions are given by

The stability of the whole calculation process is constrained by the following conditions:where Δt is the time step, ΔX = ΔY is the space step, m′ = max(Dl, m), and k′ is the correction coefficient.

The water solution of sucrose was selected as the simulation object, and the thermal property parameter values and simulation parameter values were Tm = 273.15 K, L = 333.5 J/cm3, cp = 4.21 J/cm3K, σ = 7.65 × 10−6 J/cm2, κ = 0.00131 cm2/s, Ds = 2.5 × 10−7cm2/s, Dl = 5.6 × 10−6 cm2/s, μ = 7.4 cm/Ks, ω = 2.1 × 10−4 cm, mL = 22.83, k0 = 0.075, α = 312.5, m = 0.035,  = 0.005,  = 0.1, ΔT = 0.9, θm = 0.4234, x0 = 0.56%, ΔX = ΔY= 0.005, Δt = 1.2 × 10−5, and R0 = 5ΔX.

5. Results and Discussion

5.1. Ice Crystal Growth for Interfacial Energy Anisotropy

In order to understand the influence of the strength of interfacial energy anisotropy on the growth behavior of ice crystals, Figure 6 shows the phase-field morphology of ice crystal growth corresponding to the equilibrium shape of crystals, as shown in Figure 5. When interfacial energy anisotropy is 0.04 (lower than 1/24), all orientations appear on the equilibrium shape; there is no crystal direction missing, and the main branch shape is a smooth parabola. As interfacial energy anisotropy increases to 0.06 (larger than 1/24), noise disturbance is easy to amplify, and well-developed and symmetrical side branches appear on the main stem. Due to crystal orientation disappearance, corners formed at the tips of the main stem and side branches. When is further increased to 0.3, because the influence of crystal orientation disappearance is aggravated, the main branches’ growth becomes slower and shorter and the development of side branches becomes weaker.

Figure 7 shows the relationship between the ice crystal horizontal tip temperature, solute concentration, growth rate, and radius vs. interfacial energy anisotropy. With an increase in interfacial energy anisotropy, the horizontal tip temperature first decreases and then increases and solute concentration first increases and then decreases. As a result, the growth rate first increases and then decreases and the radius first decreases and then increases. Under the dual effects of interfacial energy anisotropy promoting ice crystal growth and crystal orientation disappearance inhibiting ice crystal growth, the extreme value point of the characteristic parameters of the ice crystal tip in the horizontal direction appears at the position of an interfacial energy anisotropy strength of 0.1. It is consistent with microscopic solvability theory.

5.2. Effect of Undercooling on Ice Crystal Growth

Undercooling significantly affects the growth behavior of ice crystals. The phase-field morphology of ice crystals under different undercooling degrees is shown in Figure 8. When the undercooling degree is 0.3, the grain center extends outward to form ice crystal main branches with five-fold symmetry, and the ice crystal stem is thicker. With an increase in undercooling, the main branch becomes slender and causes necking on the root. When the undercooling degree increases to 0.7, the side branch appears on the main branch of ice crystals, which is consistent with the mechanism of crystal continuous growth in solidification theory. When undercooling continues to increase to 0.9, the main branch and side branches become more developed, forming a multibranch five-fold symmetric ice crystal structure.

Figure 9 shows the relationship between the ice crystal horizontal tip temperature, solute concentration, growth rate, and radius vs. undercooling. With an increase in the undercooling degree, the steady-state temperature at the tip of the ice crystal in the horizontal direction decreases slowly and then rapidly and solute concentration decreases rapidly and then slowly. As a result, the horizontal tip growth rate increases slowly and then quickly and the radius first decreases and then increases.

Figure 10 shows the relationship between the R2V value and undercooling. When undercooling is lower than 0.4, the value of R2V tends to be constant, and ice crystal growth is controlled by diffusion. When undercooling is larger than 0.4, the value of R2V increases gradually, and ice crystal growth is controlled by kinetics.

Finally, with an increase in undercooling, ice crystal growth changes from diffusion control to dynamic control. When undercooling reaches 0.4, the tip radius of curvature is minimized. By selecting the quick freezing conditions, undercooling is controlled at around 0.4 and the interfacial energy anisotropy at around 0.1 to inhibit the formation of large ice crystals, so as to prevent large ice crystals from damaging the cell wall, leading to the decline in the quality of quick-frozen food.

6. Summary and Conclusions

(1)Ice crystal growth morphology with five-fold symmetry can be reproduced by using the modified phase-field model. The results of the phase-field simulation are consistent with the equilibrium shape of ice crystals; the missing part of the crystal direction degenerates into a point.(2)When ɛ5 is lower than the critical value, there is no crystal direction missing, and the main branch shape is a smooth parabola. As ɛ5 increases, corners form at the tips of the main stem and side branches. Under the dual effects of interfacial energy anisotropy promoting ice crystal growth and crystal orientation disappearance inhibiting ice crystal growth, the extreme value point of the characteristic parameters of the ice crystal tip in the horizontal direction appears at the position of an interfacial energy anisotropy strength of 0.1, and the tip radius of curvature is minimized. This is consistent with microscopic solvability theory.(3)When the undercooling degree is small, the grain center extends outward to form the ice crystal main branch with five-fold symmetry, and the ice crystal trunk is thicker. With an increase of undercooling, the main branch becomes slender and the side branch appears, which is consistent with the mechanism of crystal continuous growth in solidification theory. When undercooling reaches 0.4, the tip radius of curvature is minimized, and ice crystal growth changes from diffusion control to dynamic control.(4)In this article, a phase-field model of ice crystal growth with five-fold symmetry is established, and the effects of interfacial energy anisotropy strength and undercooling on ice crystal growth morphology and tip growth behavior are discussed. During the simulation of ice crystal growth, the latent heat release of phase change was ignored, and the cell fluid was regarded as a mixed solution of water and sucrose. In the future, the influence of cooling rates, convection, and radiation heat transfer should be further considered. In view of the multiple structure of cell fluid, the ice crystal growth experiment and the simulation of the same environment were carried out by using the multiple multiphase model, which further provided a more reliable basis for improving quick freezing and fresh-keeping technology of food.

Data Availability

The data used to support the findings of this study can be obtained from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors gratefully acknowledge the financial support from the National Natural Science Foundation of China (no. 11847042), the Program of the Shangluo Science and Technology Bureau (no. SK2017-46), and the Science and Technology Innovation Team of Shangluo University (no. 18SCX004).