Numerical Simulation Study on the Front Shape and Thermal Stresses in Growing Multicrystalline Silicon Ingot: Process and Structural Design

: In this paper, a transient numerical simulation method is used to investigate the e ﬀ ects of the two furnace conﬁgurations on the thermal ﬁeld: the shape of the melt–crystal (M / C) interface and the thermal stress in the growing multicrystalline ingot. First, four di ﬀ erent power ratios (top power to side power) are investigated, and then three positions (i.e., the vertical, angled, and horizontal positions) of the insulation block are compared with the conventional setup. The power ratio simulation results show that with a descending power ratio, the M / C interface becomes ﬂatter and the thermal stress in the solidiﬁed ingot is lower. In our cases, a power ratio of 1:3–1:4 is more feasible for high-quality ingot. The block’s position simulation results indicate that the horizontal block can more e ﬀ ectively reduce the radial temperature gradient, resulting in a ﬂatter M / C interface and lower thermal stress.


Introduction
Photovoltaics (PV) is a rapidly growing market, with silicon (Si)-wafer-based PV technology accounting for approximately 95% of total production in this area in 2019. Within this field, multicrystalline (mc) technology represented approximately 34% of total production [1]. Directional solidification (DS) is an efficient method for producing mc-Si [2]. However, controlling defects represents a problem in industrial production and has restricted the development of seeded DS techniques. It is commonly known that defects are influenced by the distribution of the thermal field and the interface's front shape during the DS process [3,4]. Thermal stress is caused by nonuniform thermal deformation due to a temperature gradient. Furthermore, the front shape of a melt-crystal (M/C) interface affects crystal quality, and a convex shape is more suitable, compared with a concave shape, for excluding bubbles and impurities [5].
Since mc-Si is typically grown under sealed and high-temperature conditions, it is difficult to obtain parameters such as temperature-field distribution, crystal shape, and crystal thermal stress in the furnace through direct observation and measurement. Accordingly, simulation research has been conducted to study the effect of furnace structure and system operation on temperature distribution, front shape, thermal stress, flow field, and energy consumption. Generally, by heat and mass simulation results, adding a side insulation block in a furnace gave rise to energy saving and a smaller small-grain region [6][7][8]. However, in these simulations, the thermal stress distribution was not shown. Rao [9] showed that a narrower bottom insulation block design was preferable to one with a large horizontal width related to interface shape and lower thermal stress. Wu adjusted several hot zone parameters via simulation, including the size and position of heaters, insulation, the heat-exchange block, and the power ratio between top and side heaters to achieve an average increase in the yield rate of silicon ingots equaling 9% [10]. In addition to the insulation block, other furnace structures have also been studied. Schmid [11] used a cone-shaped crucible and susceptor for the first time to allow for lower dislocation density. Ma [12] discussed the influence of gas flow, heater position, and geometric configuration on thermal distribution and provided essential knowledge for system design. The flow pattern in the melt is another parameter that can impact silicon quality. Xie [13] showed that an enhanced flow gave rise to a more homogeneous temperature distribution among silicon melts and reduced the radial temperature gradient. Using a simulation, Stelian [14] showed that growth rate significantly influenced melt convection, and only heat flux was analyzed in the simulation. The authors of [15] showed that the melt flow pattern can be controlled by optimizing the Raleigh numbers of molten silicon during DS of the multicrystalline silicon growth process. Based on the importance of the melt pattern [16][17][18], the furnace employed a traveling magnetic field to control the convection pattern. These results indicated a higher growth rate and an improved quality of crystalline silicon ingots; nonetheless, the different effects were rather complicated. The authors of [19] showed that the melt flow and temperature distribution, particularly in the upper part of the silicon region, can be significantly influenced by a magnetic field. In these studies, the authors focused on the flow pattern in melt instead of the thermal stress in solids.
Adding insulation blocks and heat power [20] modifications are two important measures for improving the temperature field to achieve operational convenience. Although studies on the design and optimization of DS techniques have been conducted, additional efforts are needed to establish a simple modification measure and better results, and mechanism analysis should be included to offer more guidance for structure design. In this paper, a simple power ratio adjustment scheme is presented, presenting convenient implementation in current equipment without the need for modification. In addition to a horizontal block, two types of block insulation in different areas are discussed for upgrading the furnace. With simulation, the temperature distribution, front shape, and thermal stress are presented along with a detailed heat-transfer mechanism to better understand this process and serve as a reference for furnace modification. Figure 1 shows a two-dimensional axisymmetric vertical DS furnace meshed by unstructured grids. The furnace includes a quartz crucible, graphite susceptors, heat-exchange block, graphite heaters, insulation, and outside walls. A horizontal heat-exchange block was added to optimize the thermal field, leading to a low radial thermal gradient and suitable axial thermal gradient. The heaters were composed of four side heaters and one top heater. The volume of the crucible was 1100 × 550 × 550 mm 3 (length (r), width, and height (z), respectively).

Physical Model
There are two stages in the ingot growth processes: melting and solidification. In the melting process, TC1 (Figure 1) was increased to approximately 1823 K and maintained until the end of the melting process. The insulation was then slowly increased, and TC1 was decreased to approximately 1711 K to solidify. In this manuscript, we only focused on the solidification process. Table 1 lists the properties of different materials, such as density, heat capacity, thermal conductivity, and others. The silicon parameter in this table can be found in [21,22]. The crucible, plate, and insulation parameters are from the equipment company, and their emissivity was set to 0.8 for simplification.  The following assumptions were adopted to simplify the simulation processes: (1) fluid flow was ignored; (2) radiative heat exchange between solid surfaces through a nonparticipating fluid was accounted for on the assumption of gray-diffusive surface radiation; (3) heat was dispelled by cooling water in a water-cooled plate beneath the exchange block. Due to cooling, the temperature on the upper and underside surfaces of water-cooled plates remained constant and equal to 300 K, respectively. The Crysmas software package (v. 4.3.28) was used to render the simulation.
Conduction and radiation heat transfer are two primary temperature distribution types in a furnace. The temperature distribution of a DS system during the solidification process is calculated based on Fourier's fundamental laws of heat transfer. They are expressed by Equations (1) and (5), respectively.
where t is time, T is temperature, is thermal conductivity, and Q is the heat source term.  The following assumptions were adopted to simplify the simulation processes: (1) fluid flow was ignored; (2) radiative heat exchange between solid surfaces through a nonparticipating fluid was accounted for on the assumption of gray-diffusive surface radiation; (3) heat was dispelled by cooling water in a water-cooled plate beneath the exchange block. Due to cooling, the temperature on the upper and underside surfaces of water-cooled plates remained constant and equal to 300 K, respectively. The Crysmas software package (v. 4.3.28) was used to render the simulation.
Conduction and radiation heat transfer are two primary temperature distribution types in a furnace. The temperature distribution of a DS system during the solidification process is calculated based on Fourier's fundamental laws of heat transfer. They are expressed by Equations (1) and (5), respectively.
where t is time, T is temperature, k is thermal conductivity, and Q is the heat source term.
The melt-crystal interface is assumed to lie upon the melting-point isotherm, and latent heat is generated by solidification at this interface. For pure silicon with a fixed melting point, the relationship between the liquid fraction (g 1 ) and temperature is described as: where T m is the melting point of silicon.
The source term Q in Equation (1) represents the change rate of volumetric latent heat during the DS process, which can be written as where L is the latent heat of solidification.
The system primarily depended on the radiation heat transfer from the heat source to the crucible. We considered the surface to surface radiation model, calculated by Equation (5).
Heat emission from the j-th surface element: Gebhart matrix: where F is the view factors, ε i is the emissivity of the i-th surface element, and σ is the Stefan-Boltzmann constant. The basic equation for elastic stress is given below: where σ rr , σ ϕϕ , σ zz , and σ rz are the stress tensor components. An experiment was performed to verify the simulation model under the following two operating conditions. (a) No crucible was set in the furnace, and only top heating power was provided. The temperature measuring point (TC1) was above the top heater, 1 mm below the top cover, and 100 mm away from the side edge. In the simulation model, the position of TC1 was r = 540 mm and z = 1540 mm ( Figure 1). The steady and unsteady simulated temperatures of monitoring point TC1 were both compared with the experimental results. All the parameters in the simulation are the same as the experiment. For the steady case, the constant top heating power was set as 30 ± 0.6 kW, which was the same as experimental power, and the experimental and simulated temperatures at TC1 were 1705 and 1751 K (the error was −2.27%), respectively. For the unsteady case (Figure 2a), the simulated temperature showed good agreement with the experimental results when the heating power decreased from 30 to 0 kW. (b) In the operation condition, besides TC1, TC2 in Figure 1 was also monitored, and its position was r = 0 mm, z = 800 mm. In the unsteady condition, the top and side heating power both increased from 60 ± 0.2 kW to 67 ± 0.2 kW and the insulation was maintained. It can be seen from Figure 2b that, in the operation condition, the biggest temperature difference between the simulation and experiment results was 107 K (the error was −8.07%).

Effect of Power Ratio on Temperature and Thermal Stress Fields
The heat-transfer mechanisms of the top and side heaters differed from one another. The top heater transferred heat to the melt primarily based on radiation. The side heater transferred heat to the crucible's surface mainly by radiation, whereas the crucible's surface transferred heat to the melt by conduction. The heat from the top and side heaters had different effects on the radius and axial temperature gradient, which combined to determine the thermal stress, solidified front shape, and growth rate.
Four solidification stage cases were simulated to derive an efficiency ratio for the top-to-side heating power, and the total heat power equaled 52 kW ( Table 2). Case (a) (heating power ratio = 1) was set as the reference case. The furnace was cooled by cooling water (300 K) in a water-cooling jacket.  Figure 3 shows the thermal fields and thermal stresses after 16 h solidification for the four simulated cases listed in Table 2. The left side represents thermal field distribution, and the right side represents the distribution of thermal stress in ingots. Thermal stress is known as one of the major factors affecting the generation and multiplication of dislocations inside solidified silicon ingots. It is typically expressed by von Mises stress, and the distribution of thermal stress primarily depends on the temperature field. The solidification front was equal to the isothermal melting point (1685 K). By comparing temperature distribution among the four processes, the temperature fields and thermal stress distribution were clearly observed to be affected by the power ratio in all four cases, particularly in the case of temperature near the M/C interface. From Cases (a)-(d), the axial temperature gradient in the melt region decreased as top power decreased and side power increased. At the same time, the isothermal curve and the M/C interface changed from being concave to flat and even from being convex to the melt region, which was beneficial in terms of excluding impurities. In the solid region, the 1655 K isothermal line also changed from concave to flat, indicating that the axial temperature gradient in the melt-center part had decreased, whereas that in the region next to the

Effect of Power Ratio on Temperature and Thermal Stress Fields
The heat-transfer mechanisms of the top and side heaters differed from one another. The top heater transferred heat to the melt primarily based on radiation. The side heater transferred heat to the crucible's surface mainly by radiation, whereas the crucible's surface transferred heat to the melt by conduction. The heat from the top and side heaters had different effects on the radius and axial temperature gradient, which combined to determine the thermal stress, solidified front shape, and growth rate.
Four solidification stage cases were simulated to derive an efficiency ratio for the top-to-side heating power, and the total heat power equaled 52 kW ( Table 2). Case (a) (heating power ratio = 1) was set as the reference case. The furnace was cooled by cooling water (300 K) in a water-cooling jacket.  Figure 3 shows the thermal fields and thermal stresses after 16 h solidification for the four simulated cases listed in Table 2. The left side represents thermal field distribution, and the right side represents the distribution of thermal stress in ingots. Thermal stress is known as one of the major factors affecting the generation and multiplication of dislocations inside solidified silicon ingots. It is typically expressed by von Mises stress, and the distribution of thermal stress primarily depends on the temperature field. The solidification front was equal to the isothermal melting point (1685 K). By comparing temperature distribution among the four processes, the temperature fields and thermal stress distribution were clearly observed to be affected by the power ratio in all four cases, particularly in the case of temperature near the M/C interface. From Cases (a)-(d), the axial temperature gradient in the melt region decreased as top power decreased and side power increased. At the same time, the isothermal curve and the M/C interface changed from being concave to flat and even from being convex to the melt region, which was beneficial in terms of excluding impurities. In the solid region, the 1655 K isothermal line also changed from concave to flat, indicating that the axial temperature gradient in the melt-center part had decreased, whereas that in the region next to the crucible experienced little change. We thus concluded that the melt and solid center temperatures had been affected by the top heater to a larger degree compared with the side heater.

Half State
In addition, as shown in Figure 3, the largest thermal stress concentrations occurred on the corner in each case as a result of the larger temperature gradient and crucible stiffness. The maximum thermal stresses in the four cases were 3.63 × 10 8 , 3.09 × 10 8 , 1.69 × 10 8 , and 1.51 × 10 8 N/m 2 , respectively. With a decrease in the top heating power, the average thermal stress in a silicon ingot decreased, and the minimum thermal stress in the four cases decreased from 6.75 × 10 7 to 1.90 × 10 7 N/m 2 . Meanwhile, the low thermal stress region was enlarged and became more uniform. Among them, Case(a) has the largest thermal stress and uniform distribution. We thus concluded that the high top heater caused more significant thermal stress than the side heater, particularly around the M/C interface, due to the large temperature gradient of the interface.
Crystals 2020, 10, x 6 of 15 crucible experienced little change. We thus concluded that the melt and solid center temperatures had been affected by the top heater to a larger degree compared with the side heater.
In addition, as shown in Figure 3, the largest thermal stress concentrations occurred on the corner in each case as a result of the larger temperature gradient and crucible stiffness. The maximum thermal stresses in the four cases were 3.63 × 10 8 , 3.09 × 10 8 , 1.69 × 10 8 , and 1.51 × 10 8 N/m 2 , respectively. With a decrease in the top heating power, the average thermal stress in a silicon ingot decreased, and the minimum thermal stress in the four cases decreased from 6.75 × 10 7 to 1.90 × 10 7 N/m 2 . Meanwhile, the low thermal stress region was enlarged and became more uniform. Among them, Case(a) has the largest thermal stress and uniform distribution. We thus concluded that the high top heater caused more significant thermal stress than the side heater, particularly around the M/C interface, due to the large temperature gradient of the interface. This does not, however, mean that lower heat is better for the top heater. In this setup, if the top heater is decreased to 0 W, and if the solidification fraction is higher than 0.5, another interface will appear on the melt surface, which will cause a two-direction solidification. This phenomenon will occur because the temperature on the melt's top surface is too low. However, as this is not an expected condition for the DS process, it was not included in this paper.
The interface shape and temperature distribution can be clearly explained by the heat flux value and direction. Heat flux occurs in the opposite direction to the temperature gradient, whereas the temperature gradient is perpendicular to the isothermals. Because the front shape of the M/C interface was equal to the isothermal melting point (1685 K), it was affected by the heat flux direction This does not, however, mean that lower heat is better for the top heater. In this setup, if the top heater is decreased to 0 W, and if the solidification fraction is higher than 0.5, another interface will appear on the melt surface, which will cause a two-direction solidification. This phenomenon will occur because the temperature on the melt's top surface is too low. However, as this is not an expected condition for the DS process, it was not included in this paper.
The interface shape and temperature distribution can be clearly explained by the heat flux value and direction. Heat flux occurs in the opposite direction to the temperature gradient, whereas the temperature gradient is perpendicular to the isothermals. Because the front shape of the M/C interface was equal to the isothermal melting point (1685 K), it was affected by the heat flux direction at the solidification front. Figure 4 shows the heat flux distribution in the solid region when the solidification fraction increased to 50% (only half of the simulation region is shown in the figures due to the asymmetric model). In Figure 4a, the heat flux is vertical at the center before inclining toward the side in the region near the crucible. With a decrease in top power, the overall heat flux became vertical and more uniform (Figure 4). Figure 5 shows the radial and axial temperature distribution when the solidification fraction was 50%, and at that time, the interfaces were around the z = 1.0 m plane in all cases. Figure 5a shows the radial temperature distribution along the horizontal plane of z = 1.0 m. It can be seen that if the total heat was constant, larger top power led to the high temperature of the center (r = 0, T = 1690 K; Figure 5a, Case (a)), whereas a low side power led to a lower temperature on the boundary (r = 0.5, T = 1674 K; Figure 5a, Case (a)). In Figure 5a, Case (d), the temperature is as low as 1686 K at the center, whereas the temperature next to the crucible increases to 1685 K (r = 0.5). It was observed that with a top heater ratio decrease, the average radius temperature gradient of the interface decreased in the order of 0.33, 0.21, 0.07, and 0.04 K/m, respectively. Meanwhile, the axial temperature gradient also decreased. Figure 5b shows the axial temperature gradient between the horizontal plane of z = 1.0 m and the horizontal plane of z = 1.2 m. The average values of the four cases were 23.80, 11.08, 11.98, and 9.60 K/m, respectively. The maximum temperature gradient was observed in the region next to the crucible, where r = 0.5 m.
Crystals 2020, 10, x 7 of 15 at the solidification front. Figure 4 shows the heat flux distribution in the solid region when the solidification fraction increased to 50% (only half of the simulation region is shown in the figures due to the asymmetric model). In Figure 4a, the heat flux is vertical at the center before inclining toward the side in the region near the crucible. With a decrease in top power, the overall heat flux became vertical and more uniform (Figure 4). Figure 5 shows the radial and axial temperature distribution when the solidification fraction was 50%, and at that time, the interfaces were around the z = 1.0 m plane in all cases. Figure 5a shows the radial temperature distribution along the horizontal plane of z = 1.0 m. It can be seen that if the total heat was constant, larger top power led to the high temperature of the center (r = 0, T = 1690 K; Figure 5a, Case (a)), whereas a low side power led to a lower temperature on the boundary (r = 0.5, T = 1674 K; Figure 5a, Case (a)). In Figure 5a, Case (d), the temperature is as low as 1686 K at the center, whereas the temperature next to the crucible increases to 1685 K (r = 0.5). It was observed that with a top heater ratio decrease, the average radius temperature gradient of the interface decreased in the order of 0.33, 0.21, 0.07, and 0.04 K/m, respectively. Meanwhile, the axial temperature gradient also decreased. Figure 5b shows the axial temperature gradient between the horizontal plane of z = 1.0 m and the horizontal plane of z = 1.2 m.
The average values of the four cases were 23.80, 11.08, 11.98, and 9.60 K/m, respectively. The maximum temperature gradient was observed in the region next to the crucible, where r = 0.5 m.   at the solidification front. Figure 4 shows the heat flux distribution in the solid region when the solidification fraction increased to 50% (only half of the simulation region is shown in the figures due to the asymmetric model). In Figure 4a, the heat flux is vertical at the center before inclining toward the side in the region near the crucible. With a decrease in top power, the overall heat flux became vertical and more uniform (Figure 4). Figure 5 shows the radial and axial temperature distribution when the solidification fraction was 50%, and at that time, the interfaces were around the z = 1.0 m plane in all cases. Figure 5a shows the radial temperature distribution along the horizontal plane of z = 1.0 m. It can be seen that if the total heat was constant, larger top power led to the high temperature of the center (r = 0, T = 1690 K; Figure 5a, Case (a)), whereas a low side power led to a lower temperature on the boundary (r = 0.5, T = 1674 K; Figure 5a, Case (a)). In Figure 5a, Case (d), the temperature is as low as 1686 K at the center, whereas the temperature next to the crucible increases to 1685 K (r = 0.5). It was observed that with a top heater ratio decrease, the average radius temperature gradient of the interface decreased in the order of 0.33, 0.21, 0.07, and 0.04 K/m, respectively. Meanwhile, the axial temperature gradient also decreased. Figure 5b shows the axial temperature gradient between the horizontal plane of z = 1.0 m and the horizontal plane of z = 1.2 m.
The average values of the four cases were 23.80, 11.08, 11.98, and 9.60 K/m, respectively. The maximum temperature gradient was observed in the region next to the crucible, where r = 0.5 m.

End State
The end state was defined first. The melt cannot completely solidify by solely depending on the insulation lift without a power adjustment at a later period. In this case, the end state was defined as the process end state before changing power. In Case (a), the heat power had to be decreased after 10.5 h of solidification; otherwise, the solidification process would not continue. In Cases (b)-(d), the solidification time was nearly 20 h.
The temperature distribution is shown on the left in Figure 6. The maximum axial and radial temperature gradients in the solid region were (a) 4.44, 0.32, (b) 3.73, 0.14, (c) 3.42, and 0.00 K/cm. In (d), these were 3.36 and −0.2 K/cm, respectively. There was no obvious axial temperature gradient difference among these cases, whereas the temperature gradient declined as top power decreased. A negative radius temperature gradient value was observed for Case (d), meaning the interface next to the crucible was lower than that at the center. A convex interface front formed in Case (d).
Thermal stress varied throughout the entire solidification process. In the right of Figure 6, it can be observed that thermal stress reduced with a decrease in the top heat power ratio. In Cases (a) and (b), the maximum thermal stress value occurred on the periphery and at the center surface of the ingot, measured as follows: (a) 3.61 × 10 8 , (b) 2.26 × 10 8 , (c) 1.18 × 10 8 , and (d) 9.25 × 10 7 N/m 2 , respectively. Compared with those at the half stage, these values had been reduced. This indicated that among these cases, Cases (b)-(d) performed better than Case (a), and Cases (c) and (d) were more favorable for the reduction of the thermal stress value in the growth process.

End State
The end state was defined first. The melt cannot completely solidify by solely depending on the insulation lift without a power adjustment at a later period. In this case, the end state was defined as the process end state before changing power. In Case (a), the heat power had to be decreased after 10.5 h of solidification; otherwise, the solidification process would not continue. In Cases (b)-(d), the solidification time was nearly 20 h.
The temperature distribution is shown on the left in Figure 6. The maximum axial and radial temperature gradients in the solid region were (a) 4.44, 0.32, (b) 3.73, 0.14, (c) 3.42, and 0.00 K/cm. In (d), these were 3.36 and −0.2 K/cm, respectively. There was no obvious axial temperature gradient difference among these cases, whereas the temperature gradient declined as top power decreased. A negative radius temperature gradient value was observed for Case (d), meaning the interface next to the crucible was lower than that at the center. A convex interface front formed in Case (d).
Thermal stress varied throughout the entire solidification process. In the right of Figure 6, it can be observed that thermal stress reduced with a decrease in the top heat power ratio. In Cases (a) and (b), the maximum thermal stress value occurred on the periphery and at the center surface of the ingot, measured as follows: (a) 3.61 × 10 8 , (b) 2.26 × 10 8 , (c) 1.18 × 10 8 , and (d) 9.25 × 10 7 N/m 2 , respectively. Compared with those at the half stage, these values had been reduced. This indicated that among these cases, Cases (b)-(d) performed better than Case (a), and Cases (c) and (d) were more favorable for the reduction of the thermal stress value in the growth process.

Discussion
The power ratio of the top and side heaters impacted temperature distribution and silicon quality. The heater impacted the crucible's surface and silicon melt surface primarily via radiation. Equations (2)-(4) show the discrete model of total radiation flux, and the general form of radiation heat transfer between two finite surfaces, i and j, are given by Equation (10). The position of these surfaces is shown in Figure 7 (assuming that no shield existed between surfaces i and j).
where E is emissive heat, A is the area of surfaces i and j, and R is the distance between surfaces i and j.

Discussion
The power ratio of the top and side heaters impacted temperature distribution and silicon quality. The heater impacted the crucible's surface and silicon melt surface primarily via radiation. Equations (2)-(4) show the discrete model of total radiation flux, and the general form of radiation heat transfer between two finite surfaces, i and j, are given by Equation (10). The position of these surfaces is shown in Figure 7 (assuming that no shield existed between surfaces i and j).
where E is emissive heat, A is the area of surfaces i and j, and R is the distance between surfaces i and j. As shown in Figure 7, a long distance led to low heat absorption. A two-dimensional calculation was made to analyze the radiation heat exchange between the heater and the surface. Figure 8 shows the temperature distribution along the bottom horizontal plane surface when the uniform radiation heat flux was implemented on the top surface. The top heat source of 28,617 W/m 3 was used in the calculation. Figure 8 shows that even when radiation energy was equally delivered, the radius temperature distribution on the bottom surface had a strong relationship with the distance between the two surfaces. The temperature at the center of the bottom surface was higher than near the edge, and there was an approximately 80 K difference in this condition. Similarly, side-heat power caused a comparably high-temperature distribution on the crucible's surface. Based on simulation results, 1:3-1:4 is a reasonable range for the top and side power ratios, respectively.   Figure 7, a long distance led to low heat absorption. A two-dimensional calculation was made to analyze the radiation heat exchange between the heater and the surface. Figure 8 shows the temperature distribution along the bottom horizontal plane surface when the uniform radiation heat flux was implemented on the top surface. The top heat source of 28,617 W/m 3 was used in the calculation. Figure 8 shows that even when radiation energy was equally delivered, the radius temperature distribution on the bottom surface had a strong relationship with the distance between the two surfaces. The temperature at the center of the bottom surface was higher than near the edge, and there was an approximately 80 K difference in this condition. Similarly, side-heat power caused a comparably high-temperature distribution on the crucible's surface. Based on simulation results, 1:3-1:4 is a reasonable range for the top and side power ratios, respectively. Crystals 2020, 10, x 9 of 15

Discussion
The power ratio of the top and side heaters impacted temperature distribution and silicon quality. The heater impacted the crucible's surface and silicon melt surface primarily via radiation. Equations (2)-(4) show the discrete model of total radiation flux, and the general form of radiation heat transfer between two finite surfaces, i and j, are given by Equation (10). The position of these surfaces is shown in Figure 7 (assuming that no shield existed between surfaces i and j).
where E is emissive heat, A is the area of surfaces i and j, and R is the distance between surfaces i and j. As shown in Figure 7, a long distance led to low heat absorption. A two-dimensional calculation was made to analyze the radiation heat exchange between the heater and the surface. Figure 8 shows the temperature distribution along the bottom horizontal plane surface when the uniform radiation heat flux was implemented on the top surface. The top heat source of 28,617 W/m 3 was used in the calculation. Figure 8 shows that even when radiation energy was equally delivered, the radius temperature distribution on the bottom surface had a strong relationship with the distance between the two surfaces. The temperature at the center of the bottom surface was higher than near the edge, and there was an approximately 80 K difference in this condition. Similarly, side-heat power caused a comparably high-temperature distribution on the crucible's surface. Based on simulation results, 1:3-1:4 is a reasonable range for the top and side power ratios, respectively.

Furnace Structure: The Effect of the Position of the Block under the Side Heater
Based on the above results, the highest thermal stress occurred on the periphery of solid silicon. One of the measures for decreasing thermal stress is to add a block outside of the crucible to reduce the temperature gradient. When bottom block insulation is applied, the heat transfer crossing the crucible corner will decrease, which will decrease the temperature gradient. In this part, different block positions are studied to discuss the position best suited for obtaining an optimized thermal field. The positions are shown in Figure 9, i.e., (e) horizontal block, (f) a 45 • -inclined block, and (g) a vertical block. For comparison, Case (h) (without block) was also simulated. All other conditions remained the same for four cases, and the ratio of the top to the side heater was 1:3.

Furnace Structure: The Effect of the Position of the Block under the Side Heater
Based on the above results, the highest thermal stress occurred on the periphery of solid silicon. One of the measures for decreasing thermal stress is to add a block outside of the crucible to reduce the temperature gradient. When bottom block insulation is applied, the heat transfer crossing the crucible corner will decrease, which will decrease the temperature gradient. In this part, different block positions are studied to discuss the position best suited for obtaining an optimized thermal field. The positions are shown in Figure 9, i.e., (e) horizontal block, (f) a 45°-inclined block, and (g) a vertical block. For comparison, Case (h) (without block) was also simulated. All other conditions remained the same for four cases, and the ratio of the top to the side heater was 1:3.  Figure 10 shows the thermal fields (left) and thermal stresses (right) of the hot zone for three cases when the solid region reaches half the height of the crucible. Figure 10 shows that the interface gradually changes from the plane to convex to solid regions, with the block changing from horizontal to vertical. This indicates that a more appropriate S/L interface can be achieved by adding a horizontal block. Meanwhile, thermal stress gradually increased when the block changed from a horizontal to vertical position. Among four cases, the maximum thermal stress values were 1.51 × 10 8 , 1.80 × 10 8 , 3.52 × 10 8 and 4.16 × 10 8 N/m 2 , and the thermal stress distribution in Case (e) was more uniform and smaller in value compared to the other cases.

Half State
Stress distribution and the M/C interface shape depended on the temperature profiles in all cases. Figure 11 shows the radial and axial temperature distribution when the solidification fraction was 50%. Figure 11a shows the radial temperature distribution on the horizontal plane of z = 1.0 m. Case (e) exhibited a smooth temperature line and the smallest radial temperature difference (approximately 0 K) on the horizontal plane (z = 1.0 m), which provided a benefit in the form of reducing thermal stress in silicon, whereas Case (g) had the largest radial temperature difference (approximately 12 K). This result was consistent with the thermal stress simulation results. Figure  11b shows the axial temperature gradient between the horizontal plane of z = 1.0 m and the horizontal plane of z = 1.2 m, where a large temperature difference led to a high growth rate. Although Case (f) is considered as the best choice for a high growth rate, Case (e) is better suitable for high-quality products.  Figure 10 shows the thermal fields (left) and thermal stresses (right) of the hot zone for three cases when the solid region reaches half the height of the crucible. Figure 10 shows that the interface gradually changes from the plane to convex to solid regions, with the block changing from horizontal to vertical. This indicates that a more appropriate S/L interface can be achieved by adding a horizontal block. Meanwhile, thermal stress gradually increased when the block changed from a horizontal to vertical position. Among four cases, the maximum thermal stress values were 1.51 × 10 8 , 1.80 × 10 8 , 3.52 × 10 8 and 4.16 × 10 8 N/m 2 , and the thermal stress distribution in Case (e) was more uniform and smaller in value compared to the other cases.

Half State
Stress distribution and the M/C interface shape depended on the temperature profiles in all cases. Figure 11 shows the radial and axial temperature distribution when the solidification fraction was 50%. Figure 11a shows the radial temperature distribution on the horizontal plane of z = 1.0 m. Case (e) exhibited a smooth temperature line and the smallest radial temperature difference (approximately 0 K) on the horizontal plane (z = 1.0 m), which provided a benefit in the form of reducing thermal stress in silicon, whereas Case (g) had the largest radial temperature difference (approximately 12 K). This result was consistent with the thermal stress simulation results. Figure 11b shows the axial temperature gradient between the horizontal plane of z = 1.0 m and the horizontal plane of z = 1.2 m, where a large temperature difference led to a high growth rate. Although Case (f) is considered as the best choice for a high growth rate, Case (e) is better suitable for high-quality products.

End State
The same principle applies to the results shown in Figure 12. When the block position changed from horizontal to vertical, the interface front shape, isothermal line, and thermal stress behaved similarly to those in the half stage. However, the differences among them become larger, which can

End State
The same principle applies to the results shown in Figure 12. When the block position changed from horizontal to vertical, the interface front shape, isothermal line, and thermal stress behaved similarly to those in the half stage. However, the differences among them become larger, which can

End State
The same principle applies to the results shown in Figure 12. When the block position changed from horizontal to vertical, the interface front shape, isothermal line, and thermal stress behaved similarly to those in the half stage. However, the differences among them become larger, which can be represented by the maximum thermal stress values, which were 1.18 × 10 8 , 3.27 × 10 8 , 4.46 × 10 8 and 3.95 × 10 8 N/m 2 . Except for Case (e), thermal stresses in the other cases increased compared to the half state. This happened because, for the blocks in these positions, their effect on reducing heat loss became weak in the end state when the insulation setup moved higher. For this reason, the front shape became more convex when the block changed from a horizontal to a vertical position. The special front shape in Figure 12g had a significant relationship with the temperature gradient Case (g) showed a significant axial temperature gradient difference along the radius, where the smallest and largest values were 12 and 37 K/m, respectively (as shown in Figure 11b). This gave rise to the interface shape in Figure 12g, and measures will be taken to avoid this phenomenon in experiment, such as using appropriate SiN coating. It should be pointed out that Case (h) cannot solidify in this operation condition because of the low-temperature gradient after approximately 18 h. Therefore, we called this stage the end stage. This happened because, for the blocks in these positions, their effect on reducing heat loss became weak in the end state when the insulation setup moved higher. For this reason, the front shape became more convex when the block changed from a horizontal to a vertical position. The special front shape in Figure 12g had a significant relationship with the temperature gradient Case (g) showed a significant axial temperature gradient difference along the radius, where the smallest and largest values were 12 and 37 K/m, respectively (as shown in Figure 11b). This gave rise to the interface shape in Figure 12g, and measures will be taken to avoid this phenomenon in experiment, such as using appropriate SiN coating. It should be pointed out that Case (h) cannot solidify in this operation condition because of the low-temperature gradient after approximately 18 h. Therefore, we called this stage the end stage.

Discussion
To further analyze the influence of a bottom insulation block during the DS process, heat flux was also studied. The heat flux direction (in opposition to the temperature gradient and perpendicular to the isothermals), described by Equation (11), was significantly different due to the position of blocks in Cases (e)-(g).

Discussion
To further analyze the influence of a bottom insulation block during the DS process, heat flux was also studied. The heat flux direction (in opposition to the temperature gradient and perpendicular to the isothermals), described by Equation (11), was significantly different due to the position of blocks in Cases (e)-(g). q = kF(t outside − t inside ) (11) where k is the heat-transfer coefficient, W/(m·K), and F is the heat-transfer area, m 2 . Compared with the case without a block, the horizontal and inclined block (Cases (e) and (f)) changed the temperature on the crucible's surface above the block by preventing heat loss from the bottom of the side heater. Additionally, conduction heat from the crucible to the melt decreased because the temperature difference between the inside and outside of the crucible was lower, which means that (t outside − t inside ) decreased.
The vertical block (Case (g)) reduced heat from dissipating to the outside of the crucible by increasing conduction resistance around the crucible corner. Here, k is expressed by Equation (12). Adding a vertical block was equal to enlarging the thickness of the crucible (δ), which led to a smaller k value. 1 where h outside and h inside are the outside and inside of crucible surface heat-transfer coefficient, respectively (they are very small because of weak convection in this system), δ is the thickness of the crucible (0.1 m), and λ is the thermal conductivity (4.8 W/(m·K) ( Table 1)). Among three cases with blocks, when the melt solidified to 50%, the temperature difference between the melt and the outside of the crucible (the point position is shown in Figure 9) was approximately 49 (Case (e)), 53 (Case (f)), and 66 K (Case (g)), respectively. This meant that a horizontal block (Case (e)) was able to decrease heat across the crucible by approximately 28%. When the melt solidified to approximately 100%, the temperature difference between the melt and outside of the crucible was approximately 76 (Case (e)), 101 (Case (f)), and 152 K (Case (g)), respectively. This meant that Case (e) decreased heat across the crucible with the block by approximately 50%. In Case (g), approximately 20% of heat was prevented by adding a vertical block to enlarge the thickness of the crucible. Additionally, the heat decrease value was marginally influenced by the position of the insulation. Heat flux analysis, which made it easy to understand the different temperature fields in different cases, showed results consistent with the largest radial temperature gradient in the ingot.

Conclusions
Transient simulations were performed to investigate the most suitable top heat power to side-heat power ratio, as well as the effect of the block position in a modified direct solidification system. According to the results, compared with the case of a power ratio of 1:1, the most suitable temperature distribution was obtained via simulation when the power ratio was 1:3-1:4. In these cases, temperature profiles were modified and thermal stress was low and uniform, whereas the M/C interface shape changed from being convex to melt, which was beneficial for a high-quality product. Compared with the system without a block, the systems with blocks have a better performance in low thermal stress, flat front shape, low radial temperature gradient, and high axial temperature gradient. When the block position changed from vertical to horizontal, the radial temperature gradient and thermal stress decreased. Meanwhile, the M/C front shape became flatter or even slightly convex to the melt side. The simulation showed that a system with a power ratio of 1:3-1:4 and a horizontal block was particularly suitable for growing high-quality mc-Si.
Author Contributions: C.C. and G.L. contributed equally to this work as first authors. L.Z. and G.W. revised the paper. Y.H. organized figures; Y.L. contributed to the paper structure and revising sentences. All authors have read and agreed to the published version of the manuscript.