Numerical Investigation on the Effect of Thermal Gate Moving Rate on Directional Solidification Process

The temperature field distribution during the growth of crystalline silicon by the directional solidification (DS) method is an important factor affecting the growth rate, the shape of the melt-crystal (m-c) interface, and thermal stress. To solve the problem of m-c interface convexity at the early stage of crystal growth caused by supercooling at the bottom center of silicon ingot during DS. In this paper, a two dimensional (2D) global transient numerical model based on a large-size ALD-G7 (G7) crystalline silicon ingot furnace is established and experimentally verified. Based on the model, the influence of different bottom thermal gate moving process curves on the convexity of the m-c interface at the early stage was studied, with emphasis on the changes in temperature field, m-c interface, and thermal stress at the early stage of crystal growth. We have designed three cases, case 1 uses the original moving process curve of bottom thermal gate, case 2 and case 3 adjust the process curve to 0.95 and 0.9 of the original ratio, respectively. The numerical simulation results show that compared with case 1, the average heat dissipation at the bottom of the silicon ingot in case 2 and case 3 is reduced by 2.7 kW·m−2 and 4 kW·m−2, and the maximum thermal stress is also reduced by 3 MPa and 4 MPa, separately. The maximum convexity of the interface is reduced by 25 mm and 20 mm, and the convexity is reduced by 55% and 44% on average, respectively.


Introduction
Polycrystalline silicon has dominated the photovoltaic (PV) market for more than 10 years due to its economic and electrical performance [1]. However, according to the ITRPV report [2], monocrystalline silicon will dominate the market in the future. Directional solidification growth of crystalline silicon has faced unprecedented challenges. To enhance the market competitiveness of the directional solidification method for silicon growth, and improve its electrical properties, production efficiency, and photoelectric conversion efficiency, it has become the main direction of future research on directional solidification method [3][4][5]. Impurity, dislocation, random grain boundary, and other defects seriously affect the quality of silicon ingot. How can we control the defects has always been a difficult problem in industrial production [6,7]. These defects can be controlled by optimizing the m-c interface and thermal stress, so optimizing the temperature field and process parameters in the furnace from a macroscopic point of view is very important to improve the environment of crystal growth [8].
To optimize the m-c interface and reduce thermal stress, Sundaramahalingam et al. designed a crucible with a concave bottom to strengthen the local cooling at the bottom of the crucible and obtained a better m-c interface by comparing the four thicknesses of the bottom [32]. Yang et al. added a conical heat insulation device at the bottom of the traditional furnace to obtain crystals with lower convexity, lower thermal stress, and better quality compared with the traditional furnace [33]. Keerthivasan et al. compared and analyzed the crystal growth process of conventional retort and retorts with insulation materials inserted at the top, middle, and bottom, respectively [34]. Then they found that the m-c interface of the crystals grown with insulation material inserted in the middle of the retort was slightly convex, with low thermal stress, and received the highest quality. Su et al. designed pyramid bottom grille based on the original bottom grille, which improved the supercooling in the central region and the m-c interface at the early stage of crystal growth, and reduced the thermal stress in the bottom center and corner by 20.22 and 16.14%, respectively [35].
The studies above have proved that the reasonable temperature field and process design did improve the temperature field of crystal silicon growth, so that crystal growth can obtain a more uniform temperature field and less thermal stress, and the realization of a slightly convex m-c interface is conducive to improving the quality of silicon ingot [21,36]. To better control the temperature field and effectively affect the shape of the m-c interface, so as to improve the quality and yield of silicon ingot, the shape of the m-c interface at the early stage of crystal growth has a profound influence on the growth of later stage. And since the m-c interface at the early stage of crystal growth is affected by the heat dissipation at the bottom, so we decided to pay attention to the moving rate of the bottom thermal gate. Due to the temperature field is difficult to obtain through experiments, this paper adopts the method of numerical simulation to propose three cases with different bottom thermal gate moving rates, compares and discusses the temperature field, flow field, m-c interface, flow velocity, heat and flow distribution, thermal stress distribution, oxygen (O) and carbon (C) impurities in different cases at the early stage of crystal growth.

DS System Description
For crystal silicon growth, three dimensional (3D) numerical simulation can predict temperature field, flow field and m-c interface more accurately [9]. However, due to the high computational cost and difficulty in the implementation of the 3D model, a simplified 2D axisymmetric model was adopted as confirmed by many researchers [15,21,37,38]. The 2D model we used in this paper was verified by comparing simulated and experimental temperature and power changes. Please refer to our previous research for details [35,39]. The G7 furnace model used in numerical simulation calculation is shown in Fig. 1, including top, side, and bottom insulation, top, side, and bottom heater, cover, quartz crucible, susceptor, heat exchange block, grille, thermal gate, and cold copper plate. The size of the silicon ingot is 1320 × 1320 × 385 mm 3 . After the silicon material completely melted, the heat dissipation mainly depends on the grille, cold copper plate, and thermal gate under the bottom heater.
The detailed description of heat and mass transfer models and the governing equations of heat transfer, flow, and melt convection used in the model, the boundary conditions and detailed transport equations of transportation of O and C impurities, material property parameters, and its simplifications can be found in our previous studies [39][40][41].
The analysis of thermal stresses in silicon ingot is applied an axisymmetric displacement-based thermo-elastic stress model. The governing partial differential equations for momentum balance in axisymmetric models can be described as follows [42]: The stress-strain relationship of the crystal is taken as: where σ ij , C ij , ε ij , E, γ are stress tensor, elastic matrix, strain tensor, Young's modulus, Poisson's ratio. β = 4.5 × 10 −6 K −1 , E = 165 GPa, γ = 0.29. There is no traction boundary � ⃗ • � ⃗ n = 0 . Thermal stress is usually expressed by the Von Mises stress. The Von Mises stress σ von is defined as: The temperature of the shell of furnace is set as the normal temperature boundary of 300 K; The flow rate of argon gas inlet is determined according to the furnace type and the operation parameters of growth stage. Argon gas outlet adopts natural flow boundary condition.
The number of grids used for the model calculation was 32,580, and the grids were divided by ICEM, using a mixture of structured and unstructured grids. Regular areas such as silicon area, quartz crucible, susceptor and heat exchange block are structured grids, and irregular areas such as argon area are unstructured grids. The open and close of the bottom thermal gate involved in this paper are realized by dynamic mesh. There are three mesh updating methods in Fluent, namely smoothing, layering and remeshing, which are used to move boundary and complete mesh updating to avoid negative volume grids.

Bottom Thermal Gate Moving Rate Design
Based on the previous analysis of the heat zone parameters in G7 furnace, we found that there is a problem of supercooling at the bottom center of crystalline silicon growth in the furnace, which leads to the over-convex m-c interface at the early stage of crystal growth. The over-convex m-c interface at the early stage will inevitably affect the m-c interface at the later stage. The growth rate of the m-c interface at the early stage is mainly related to the moving rate of thermal gate. The faster the moving rate, the faster the heat dissipation, the faster the crystal growth of the center, which eventually leads to the problem of over-convex m-c interface at the early stage. In order to study the effect of the change of the G7 bottom thermal gate moving process curve on the growth of silicon at the early stage of crystal growth, we decided to adjust the bottom thermal gate moving rate in the first 8 hours of crystal growth. Figure 2 shows the bottom thermal movement design. Among these, case 1 is the original unadjusted bottom thermal gate process curve, case 2 and 3 adjusted the bottom thermal moving distance to 0.95 and 0.9 under the original ratio, respectively. Figure 3 shows the temperature distribution and flow field structure under different bottom thermal gate movement rates at the 8 hours stage of crystal growth process. At the bottom of the center of the silicon ingot, that is, just above the opening of thermal gate, the bottom temperature of the silicon ingot for case 1 to case 3 are 1621 K, 1623 K, and 1635 K, respectively. Since the open direction of the bottom thermal gate is from the center to two sides, the bottom of the silicon ingot, especially the center, will be cooling first. And at the same time, the extent of cooling here is directly related to the moving rate of the thermal gate. According to our design, the rate of thermal gate: case 1 > case 2 > case 3, so the same 8 hours at the time, the opening width of bottom thermal gate: case 1 > case 2 > case 3. Therefore, at this time, the bottom from case 1 is the most severe cold, which also results in the lowest temperature at the bottom. Figure 4 shows the temperature distribution of the side wall of the silicon ingot along the height of silicon ingot for Fig. 2 The design of bottom thermal gate for moving distance Fig. 3 The temperature field (left) and flow structure (right) of the silicon area at 8 hour of crystal growth process for (a) case 1; (b) case 2; (c) case 3 the three cases. As can be seen from the figure, the overall temperature of case 3 is higher than that of case 2 and case 1, followed by case 2, and the overall temperature of case 1 is the lowest, showing a decreasing trend from case 3 to case 1. The average temperature of the ingot side wall from case 1 to case 3 are 1699 K, 1701 K, and 1704 K, respectively. Figure 5 shows the temperature distribution along the radial distance of the bottom wall of silicon ingot for the three cases. The temperature distribution of the three cases in Fig. 5 is basically consistent with that in Fig. 4. The average temperature of the silicon bottom wall from case 1 to case 3 are 1646 K, 1650 K, and 1657 K, respectively. Compared with the side wall, the position of the bottom wall is much closer to the thermal gate, so after the thermal gate opens for heat dissipation, the temperature of the bottom wall is much lower than the side wall. From the point of view of the average temperature difference, the maximum average temperature difference of the side wall is only 5 K, while the average temperature difference of the bottom wall is 11 K. The average temperature difference of the bottom side wall from case 1 to case 3 are 53 K, 51 K, and 47 K, respectively. As the lateral movement rate of bottom thermal gate decreases and the heat dissipation rate slows down, the temperature at the side and bottom of the silicon ingot are higher than that in case 2 and case 3.

Comparison of Temperature and Flow Filed
As shown in Fig. 3, the distribution of flow field structure in the three cases are similar. They all have two vortex cells: A and B, with the same flow direction and counter-clockwise flow. The maximum flow rate of melt at the m-c interface in three cases: case 1 > case 3 > case 2. In the process of crystal growth by directional solidification method, the melt flow is driven by the radial temperature difference of the melt. The radial temperature difference for the three cases at 10 mm above the m-c interface is 2.2 K, 1.52 K, and 1.98 K, respectively. From the distribution of flow field at the m-c interface shown in Fig. 6, it can be found that the comparison of the maximum flow velocity for three cases is consistent with that shown in Fig. 3: case 1 > case 3 > case 2, but the flow rate at the side wall of silicon ingot presents as case 3 > case 2 > case 1. The maximum velocity for case 1 is the largest, but the average velocity of the three cases is 2.1E-3 mm·s −1 , 1.7E-3 mm·s -1, and 2.2E-3 mm·s −1 , respectively. Figure 7 shows the comparison diagram of m-c interface shape of crystal growth at 8 hours for the three cases. From the figure, we can intuitively see that the convexity of the interface for the three cases is: case 1 > case 3 > case 2. Compared with case 1 with the maximum convexity of interface, case 2 and case 3 reduced the convexity by 55% and 44% on average, respectively. From the perspective of cooling, with the slowing down of the moving rate of the bottom thermal gate, the cooling situation at the center of the silicon ingot is improved very well, and the middle convex situation will not be caused by the fast heat dissipation in the middle of silicon ingot or the fast crystal growth.

Comparison of m-c Interface
The moving rate of the thermal gate is also reduced and the rate for case 3 is slightly lower than that of case 2, but the convexity of case 3 is greater than that of case 2. There are two reasons for the convexity of the m-c interface. One is that the silicon ingot grows too fast in the middle, and the other is that the silicon ingot grows too slowly on both sides. Although case 3 has the slowest moving rate of the thermal gate, the central crystal growth of silicon ingot should also be the slowest, but we can not ignore the influence of side heat flow on the side crystal growth rate.

Comparison of Heat Flux Profiles
In order to better explore the difference in side crystal growth, it is necessary to reflect the change of heat flux on the side wall more intuitively. Figure 8 shows the distribution of heat flux size on the side wall of silicon ingot. It can be found as follows: case 3 > case 2 > case 1. The average heat flux value of the side wall from case 1 to case 3 is 1.5 kW·m −2 、2.07 kW·m −2 and 2.4 kW·m −2 , respectively. Compared with case 1, case 2 and case 3 increased by 38% and 60%. At the m-c interface, the heat flux value of case 2 and case 3 is positive, indicating that there is no edge warping at the side of m-c interface at this time, while heat flux for case 1 is negative and there is edge warping. Figure 9 shows the heat flux distribution of the bottom wall of silicon ingot. The heat dissipation of the bottom wall of silicon ingot is as follows: case 1 > case 2 > case 3. The average heat dissipation from case 1 to case 3 is 14.6 kW·m −2 、11.9 kW·m −2 Fig. 7 Comparison of meltcrystal interface shape at 8 hour for crystal growth process Fig. 8 Heat flux profiles along the silicon side wall at 8 hour for crystal growth Fig. 9 Heat flux profiles along the silicon bottom wall at 8 hour for crystal growth and 10.6 kW·m −2 , respectively. The heat dissipation of case 2 and case 3 is 18% and 27% lower than that of case 1.
The changes of the heat flux are caused by the moving rate of the thermal gate. With the decrease of the moving rate, the opening size of gate at the same time decreases, and the loss of heat flux decreases. The temperature of the side wall will continue to remain high if the heat dissipation from the bottom is not effective, resulting in high heat flux of the side wall. It is worth noting that the heat flux at the bottom center of case 3 is slightly higher than that of case 2, which we believe this may be caused by the the flow velocity affected by the small opening size. It also reminds us that simply slowing down the moving rate of the thermal gate is not an absolute beneficial approach. From case 1 to case 3, the heat dissipation from the bottom of the silicon ingot is gradually reduced, and the growth rate of the silicon ingot center is slowing down, which plays a role in restraining the central convexity of the m-c interface, and effectively improves the convexity of the m-c interface. Figure 10 shows the temperature field and thermal stress distribution of silicon crystal in three cases. The results show that there are large thermal stress values in the bottom center and bottom corner of the silicon ingot in all three cases, which means that the radial temperature gradient at the bottom is larger than that at the top. The maximum thermal stress from case 1 to case 3 is 15 MPa, 12 MPa, and 11 MPa, respectively. Compared with case 1, the maximum thermal stress of case 2 and case 3 is reduced by 20% and 27%, respectively. In addition, the minimum thermal stress of the three cases is 3.2 MPa, 2.2 MPa, and 1.8 MPa, respectively. The minimum thermal stress of case 2 and case 3 is reduced by 31% and 43%, respectively. The gate opening rate decreases and the open width decreases, results in the decrease of cooling rate at the bottom of the crucible, and the temperature gradient decreases as the temperature difference decreases, which effectively reduced the thermal stress in the silicon ingot. Figure 11 shows the thermal stress distribution along the radial distance at the bottom of the silicon ingot. The peak value appears at the center and the corner of the silicon ingot, and the thermal stress at the bottom of the silicon ingot from case 1 to case 3 also shows a descending trend.

Comparison of Impurities Distributions
As can be seen from Fig. 12, the overall distribution form of O and C impurities in the silicon region from case 1 to case 3 is not very different. The concentration of O impurities near the side wall of the crucible is relatively high, and the concentration at the center of the melt and the free surface of the melt is relatively low. C impurity mainly concentrates at the center of the free surface and near the m-c interface. The former is because it is located in the dead water zone where silicon melt flows, and the latter is because of the separation effect at the m-c interface.
The concentration of C impurity near the m-c interface is as follows: case 3 > case 1 > case 2, it is affected by the flow velocity near the m-c interface. The C impurity near the free surface of the melt in case 3 is significantly larger than that in case 1 and case 2. Since the thermal gate at the bottom  moves more slowly in case 3, the sidewall of crucible lasts longer under high temperature and the heat dissipates more slowly. Under the condition of high temperature, more oxygen impurities produced from the sidewall at the same time, then the increase of oxygen impurities evaporated from the m-c interface will lead to more CO flow back to the interface, bring more C impurities into the melt. And the flow velocity of vortex cell A on the free surface from case 3 is weaker than that in case 1 and case 2 (see Fig. 3), the transportation of impurities is suppressed, so the concentration of C impurities near the free surface from case 3 is significantly higher than that of case 1 and case 2.

Conclusions
In this paper, the temperature field, flow field, thermal stress field, and impurities in the silicon region of crystal growth for different bottom thermal gate moving rates in G7 furnace are analyzed. With the decrease of moving distance of thermal gate from case 1 to case 3, the flow structure and strength changed little, but the cooling condition of the center of silicon ingot gradually improved, and the convexity of interface is reduced by 55% and 44% on average, respectively. From case 1 to case 3, due to the gradual reduction of the axial temperature gradient, the maximum thermal stress of case 2 and case 3 is reduced by 20% and 27%, respectively.
In conclusion, the cases we proposed to optimize the moving process curve of bottom thermal gate can achieve the purpose of improving the over-convex m-c interface at the early stage of crystal growth and reducing the thermal stress of silicon ingot.