Simulation Analysis of Silicon Ingot Growth in Directional Solidification System

1College of Artificial Intelligence, Yango University, Mawei District, Fujian 350015, China 2Department of Mechanical Engineering, National Kaohsiung University of Applied Science, Kaohsiung 807, Taiwan 3Department of Chemical and Materials Engineering, National University of Kaohsiung, Kaohsiung 811, Taiwan 4Department of Aeronautical Engineering, Chaoyang University of Technology, Taichung 413, Taiwan


Introduction
Recently, more domestic and foreign researchers and companies have focused their research on the growth of polycrystalline silicon. The investigation has focused on how to reduce the impurity concentration and enhance the growth rate and on analyzing the distribution of the thermal field in the growth furnace (1)(2)(3)(4)(5) with the aim of growing silicon ingots with higher solar photovoltaic conversion efficiency that can be used to fabricate solar cells with greater costeffectiveness. A directional solidification system (DSS) is mainly used to manufacture and produce polycrystalline silicon, and most manufacturers use the DSS furnace produced by GT Solar Ltd. The DSS furnace allows axial heat conduction and minimizes heat loss in the horizontal direction, enabling the silicon ingot to grow in one direction. Because the one-round growth of a polycrystalline silicon ingot takes a long time, making its experimental investigation time-consuming, we have developed a numerical analysis method to simulate and analyze the distribution of the thermal field in the growth furnace. We improved the internal structure of a furnace and changed the state of the temperature field to obtain the optimal growth parameters. (4)(5)(6)(7)(8) We analyzed the internal temperature field distribution and discussed the ingot growth process in the DSS furnace. We compared the results obtained from numerical analysis and an actual experiment and improved the growth process. When the simulation technology is used to find the optimum parameters, less labor, materials, and time are required for the manufacture of silicon ingots, and the simulation parameters of the ingot growth furnace can be changed arbitrarily without adding to the cost.
Numerical analysis technologies are widely used to simulate the cooling rate of silicon ingots, and computational fluid dynamics (CFD) is now one of the most widely used technologies in the analysis of fluid mechanics. The finite volume method (FVM) is a discrete technology that has rapidly developed in recent years. The FVM has high computational efficiency and is widely used in the CFD field. Many CFD software packages have been developed to use the FVM, for example, Fluent, ANSYS CFX, and STAR-CD. Fluent is one of the most widely used CFD software packages because it can simulate and analyze the fluid flow and heat transfer in complicated geometric areas. Fluent allows the use of flexible mesh characteristics by supporting many different mesh structures; thus, users can freely choose unstructured or structured meshes to divide the geometric area. For example, triangle and quadrilateral meshes can be used to simulate two-dimensional structures, and tetrahedron, hexahedron, cuneiform, orthorhombic pyramid, and polyhedron meshes can be used to simulate three-dimensional structures. Fluent can also support the use of mixed meshes and has a user-friendly interface that is easy for new users to use. Fluent was designed on the basis of the thinking of a CFD software group, i.e., from the perspective of users. Using Fluent, we can focus on the physical phenomena of complex fluid fields, because this software has the functions of different discrete grids and specialist numerical methods to obtain the best combination of calculation speed, stability, and accuracy. Fluent can be used to solve complex fluid calculation problems in different research fields with high efficiency.
Silicon ingots are used as semiconductors in many applications and fields, including the chemical and physical sensors used in genetic diagnostics, pollution control, medicine, biophysics, biochemistry, asthma diagnosis, drug delivery, optical switching, and optical sensors. (9,10) Therefore, the manufacture of high-quality ingots for the substrates of sensors is very important. In this study, we examined silicon ingot growth in a DSS furnace. (6,(11)(12)(13)(14) Next, we used a numerical analysis method to find the effects of different parameters on the solidliquid phase distribution map in the silicon ingot growth process in a DSS furnace, including the geometric model, the meshes, the governing equations, and the calculation processes and boundary conditions of the used algorithm.

Simulation Parameters
The simulation analysis of silicon ingot growth in a DSS furnace using software was divided into three processes: pretreatment, numerical calculation, and postprocessing.
1. Pretreatment: we used Autodesk computer-aided design (AutoCAD) software to construct and lay out the two-dimensional model and used the design modeler to repair the model, specify boundaries, and finally grid the meshes. 2. Numerical calculation: we used Fluent to perform the numerical calculation, along with a visual studio functional library to designate the user's defined functions. 3. Postprocessing: we imported the results calculated using Fluent to perform the data processing of the tabular ribbon distribution map, contour distribution map, and point data.
The silicon ingot growth involves five processes: heating, melting, growth, annealing, and cooling. Figure 1(a) shows the two-dimensional layout of the main parts and Fig. 1(b) shows the work flow field of the silicon ingot growth furnace manufactured by GT Solar Ltd.; we used this DSS furnace for the simulation of heating fields. The temperature in the ingot growth furnace was controlled by the heater, and there were two temperature sensors, TC1 and TC2. TC1 was used to monitor the temperature of the heater, and TC2 was used to monitor the temperature of the graphite protection cage. A water-cooling system was included in the furnace wall, and the temperature outside the furnace was kept at room temperature. We used AutoCAD to construct a two-dimensional global-model mesh, as shown in Fig. 1 Skewness was used to evaluate the mesh quality and is defined as where G is the optimal cell size and H is the cell size. Skewness has a value between 0 and 1.
When the skewness is 0, the mesh has the optimal characteristic and consists of regular polygons (equilateral triangles or squares). A skewness of 1 suggests that the differences between the maximum and minimum angles of the polygon meshes is too large. Table 1 shows the skewness values for different qualities of the mesh. The default setting of Fluent is that when the skewness value is larger than 0.97, the calculation process is terminated. However, the skewness of the initial meshes and remeshes must be controlled to less than 0.97, and the smaller the skewness, the better. The insulation cage was divided into movable and unmovable regions, as shown in Fig. 2. The unmovable region used triangle unstructured meshes, allowing the region to go through the remeshing process. There are two remeshing mechanisms: one evaluates the length of the meshes and the other evaluates the skewness. When the length of one mesh was smaller than 6.5 mm, it was combined with neighboring meshes, and when the length of one mesh was larger than 11 mm, it was split into two meshes. When the skewness of one mesh was larger than 0.6, it was combined with neighboring meshes or split to reduce the skewness to less than 0.6. The other parts and the fluid field were constructed using quadrilateral structured meshes. Figure 3 shows schematic diagrams of the initial mesh and the remeshes after different times.
We first compared the simulated numerical calculation results of the DSS ingot furnace with the actual results of the fabrication process. The simulated results were compared with the temperature extracted from the TC2 temperature sensor. Also, the temperature changes were

Simulation Processes and Results
In this study, the FVM was used to find the solution of the equations governing the DSS furnace. The governing equations of each control volume are as follows. First, the continuity equation is where t is time, ρ is density, u x is the velocity component in the x direction, u y is the velocity component in the y direction, and u z is the velocity component in the z direction. In an uncompressed fluid, ρ is constant and Eq. (2) can be rewritten as When the Hamiltonian differential operator is introduced into Eq. (3), the equation becomes which can also be expressed as Next, the momentum equations are as follows: where P is the static pressure, τ is the component of shear stress, and f is the mass force. If the momentum is only affected by the gravity force and the y-axis is vertically downward, then f x = f z = 0 and f y = g. The energy equations can be expressed as [ ] where ( )  , and S depend on the time, convection, conduction, material diffusion, viscous dissipation, and enthalpy. E is the total energy (internal energy + dynamic energy + static energy), u  is the velocity vector, h j is the enthalpy, J j is the diffusion of the fluid, k eff is the thermal conductivity, and S is the heat source item.
Next, the turbulence equations (k-ε equations) are given by i n wh ich 1 max 0.43, where a is the gas absorption coefficient, σ s is the gas scattering coefficient, I is the light intensity, r  is the position vector, s  is the direction vector, s′  is the scattering direction vector, T is the temperature of the local gas, σ is the Boltzmann constant (5.669 × 10 −8 /m 2 ·K 7 ), Φ is the phase function, ′ Ω is the solid angle, and n is the wavelength. Finally, the energy equations of the solidification problem are � and where H is the enthalpy, v  is the velocity vector of flow field, ΔH is the latent heat, c p is the specific heat, h ref is the reference enthalpy, and T ref is the reference temperature.
To show the difference between the simulation and actual conditions, Fig. 5 shows the measured and simulated temperature curves at point TC2. The error tolerance is about 6%. This result demonstrates that the CFD software can be used to simulate the process of silicon ingot growth in the DSS furnace. Next, we used the variation of the slope to find the change in the solid-liquid interface shape. We used the lowest point as a benchmark, the section from the lowest point to the crucible-wall area was referred to as the wall-face section, and the highest point of the crucible-wall section was located on the crucible wall. Using the lowest and highest points, we calculated the average slope in the solid-liquid interface of the crucible-wall section, which we called the concavity, the section from the lowest point to the center of the crucible was referred to as the central section, and in the central section, the highest point of the solid-liquid interface was usually at the center of the molten silicon. The average slope of the solid-liquid central section calculated from the highest and lowest points was referred to as the convexity. The concavity was defined as a negative value, and the smaller its value, the more concave the  Fig. 4, where the vertical spacing between each line is the temperature difference between different points. A larger spacing indicates that the vertical temperature gradient is higher, meaning that the DSS furnace is more able to drive the growth of the crystal silicon vertically upwards. Among points A, B, and C, point C, at the highest location and nearest the heater, had the lowest cooling rate, and 36 h was required to completely release the latent heat. Point A had the lowest location and the highest cooling rate, and only 0.6 h was required to completely release the latent heat. Therefore, the growth rate of the silicon ingot was highest at the beginning, then slowly decreased. Points A*, B*, and C* were closer to the crucible-wall section than points A, B, and C, but their cooling rates were higher than those of points A, B, and C in the central section of the molten silicon. Figure 6(a) shows the temperature difference curve between points A and C, which is the axial temperature change. When the silicon at point A solidified, point C had a high temperature of 1685 K, and the silicon at this point was still molten. The heat transfer characteristic of the solidified silicon was better than that of the molten silicon, and the silicon at the bottom of the furnace had the fastest heat dissipation. The temperature difference between points A and C increased gradually until the silicon at point C was solidified. As shown in Fig. 6(a), the maximum temperature difference between points A and C was 88.41 K. Because the silicon at point C was solidified after 36 h, the temperature decreased gradually.  Figure 6(b) shows the difference in the radial temperature between each point in the central section and the corresponding point in the crucible-wall section. The change in the radial temperature in the top area was because the Ar gas flow had a strong influence on temperature gradient, and the silicon at point C was solidified first. Point C* was closer to the heater than point C, and its temperature remained at 1685 K; thus, the temperature difference between the two points increased with time. The difference in the central radial temperature resulted in the silicon at point B* being solidified at 7 h and the silicon at point B being solidified at 14 h, and the maximum temperature difference between the two points was 17.84 K. When the silicon at point B was solidified, the difference in the radial temperature dropped. The difference in the central radial temperature at the bottom meant that the silicon at point A* was instantly solidified at the beginning of growth and the silicon at point A was solidified at 0.6 h, then the difference in the radial temperature dropped gradually.
Diagrams showing the variation of the temperature profile as a function of cooling time are shown in Fig. 7. As the cooling time increased from 10 to 40 h, it can be clearly seen that the bottom of the silicon gradually cooled. The thermal energy escaped from the bottom clearance of the graphite insulating plate, and the radial temperature gradient of the silicon material in the crucible was proportional to the distance from the thermal heater. Figure 8 shows the simulated solid-liquid phase distribution map, where the cyan zone was the region of 100% solid silicon, the red zone was the region of 100% liquid silicon, and the zone between the cyan and red zones was the solid-liquid fuzzy zone. To extract the required points with the lowest concavity and convexity, the area of the 50% fuzzy zone was defined as the solid-liquid interface. Figure 9 shows the solid-liquid phase distribution maps at different times. Because the crucible side walls cooled fastest, the silicon ingot started to grow from the two side walls. If the radial temperature gradient can be well controlled to ensure that it is not too large, we can prevent the silicon ingot from starting to grow from the side walls. If the silicon ingot grows from the side walls, a long-strip-type crystal lattice will form. Such a crystal lattice would cause the dopants to accumulate in the grain boundaries, and the higher impurity concentration would degenerate the energy conversion efficiency of the grown silicon wafer.
As shown in Fig. 9, as the cooling time increased from 5 to 15 h, the solid-liquid interface changed from flat to convex, the convexity increased, and the solid silicon had a convex outward appearance. With increasing cooling time, the crystal grains of the central section grew more outward, and the growth direction of the silicon ingot was not sufficiently vertical. As shown in Fig. 10, the maximum convexity was 0.23 and the minimum concavity was −0.02. These simulation results have demonstrated that we can use CFD software to simulate the growth conditions of a silicon ingot in a DSS furnace. We believe that simulation technology can optimize the growth parameters of the directional solidification process for growing polycrystalline silicon ingots with higher efficiency.

Conclusions
In our simulation of crystal growth in a DSS furnace, we used triangle unstructured meshes for the unmovable region and quadrilateral structured meshes for the other parts and the fluid field. The continuity equation, the momentum equations, the energy equations, the turbulence equations, and the radiation equation were successfully used to simulate the growth of a silicon ingot in a DSS furnace. We found that the maximum temperature difference between points B and B* was 17.84 K. From the simulation results, we also found that as the cooling time increased from 5 to 35 h, the maximum convexity was 0.23 and the minimum concavity was −0.02. These results demonstrated that the simulation technology can optimize the growth parameters of the directional solidification process for growing polycrystalline silicon ingots with higher efficiency.