Next Article in Journal
Influence of Al and N Content and Cooling Rate on the Characteristics of Complex MnS Inclusions in AHSS
Next Article in Special Issue
Pyroelectric Properties of Bismuth Borate BZBO Single Crystals
Previous Article in Journal
A Coloring Study of the Ga Richest Alkali Gallides: New In- and Hg-Containing Gallides with the RbGa7- and the K3Ga13-Type Structure
Previous Article in Special Issue
Spectroscopic, Raman, EMPA, Micro-XRF and Micro-XANES Analyses of Sulphur Concentration and Oxidation State of Natural Apatite Crystals
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Energy Institute, Qilu University of Technology (Shandong Academy of Sciences), Jinan 250100, China
2
School of Energy and Power Engineering, Qilu University of Technology (Shandong Academy of Sciences), Jinan 250100, China
3
State Key Laboratory of Crystal Materials, Shandong University, Jinan 250100, China
*
Author to whom correspondence should be addressed.
Chengmin Chen and Guangxia Liu contributed equally to this work as first authors.
Crystals 2020, 10(11), 1053; https://doi.org/10.3390/cryst10111053
Submission received: 13 October 2020 / Revised: 12 November 2020 / Accepted: 15 November 2020 / Published: 19 November 2020

Abstract

:
In this paper, a transient numerical simulation method is used to investigate the effects of the two furnace configurations on the thermal field: the shape of the melt–crystal (M/C) interface and the thermal stress in the growing multicrystalline ingot. First, four different 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 flatter and the thermal stress in the solidified 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 effectively reduce the radial temperature gradient, resulting in a flatter M/C interface and lower thermal stress.

Graphical Abstract

1. 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.

2. Physical Model

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 mm3 (length (r), width, and height (z), respectively).
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.
ρ C p T t = · ( k T ) + Q
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 (g1) and temperature is described as:
g l = { 1 , when   T > T m 0 , when   T < T m
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
Q = ρ ( Δ H ) t
Δ H = g l L
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).
q i i n = j = 1 N G j i q j o u t
Heat emission from the j-th surface element:
q j o u t = ε j σ T 4
Gebhart matrix:
G = I ε F ( I ( 1 ε ) F ) 1
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:
1 r ( r σ r r ) r + σ r z z σ φ φ r = 0
1 r ( r σ r z ) r + σ z z z = 0
where σ r r , σ φ φ , σ z z ,   and   σ r z 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%).

3. 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.

3.1. Half State

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.
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 × 108, 3.09 × 108, 1.69 × 108, and 1.51 × 108 N/m2, 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 × 107 to 1.90 × 107 N/m2. 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 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.

3.2. 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 × 108, (b) 2.26 × 108, (c) 1.18 × 108, and (d) 9.25 × 107 N/m2, 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.

3.3. 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).
Q = ( E b 1 E b 2 ) c o s φ i · c o s φ j π R 2 δ A j δ A x
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/m3 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.

4. 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.

4.1. Half State

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 × 108, 1.80 × 108, 3.52 × 108 and 4.16 × 108 N/m2, and the thermal stress distribution in Case (e) (Figure 10a) was more uniform and smaller in value compared to the other cases.
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.

4.2. 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 × 108, 3.27 × 108, 4.46 × 108 and 3.95 × 108 N/m2. Except for Case (e) (Figure 12a), 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 12c 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 12c, 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) (Figure 12d) cannot solidify in this operation condition because of the low-temperature gradient after approximately 18 h. Therefore, we called this stage the end stage.

4.3. 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 = k F ( t o u t s i d e t i n s i d e )
where k is the heat-transfer coefficient, W/(m∙K), and F is the heat-transfer area, m2.
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 (toutsidetinside) 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 k = 1 h outside + δ λ + 1 h inside
where houtside and hinside 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.

5. 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.

Funding

This research was funded by Shandong Province Natural Science Foundation (ZR2018BEM027, ZR2017LEE018), Shandong Double Hundred Talent Plan (WST2017001), Open Fund of State Key Laboratory of Crystal Materials (Shandong University, No. KF1911), the Pilot Project Scheme of Shandong Academy of Sciences, and Education and Industry in Qilu University of Technology (Shandong Academy of Sciences) (2020KJC-YJ04).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Fraunhofer Institute for Solar Energy Systems. Photovoltaics Report; Fraunhofer ISE: Freiburg, Germany, 2020. [Google Scholar]
  2. Li, T.F.; Huang, H.C.; Tsai, H.W.; Lan, A.; Chuck, C.; Lan, C.W. An enhanced cooling design in directional solidification for high quality multi-crystalline solar silicon. J. Cryst. Growth 2012, 340, 202–208. [Google Scholar] [CrossRef]
  3. Chen, X.J.; Nakano, S.; Liu, L.J.; Kakimoto, K. Study on thermal stress in a silicon ingot during a unidirectional solidification process. J. Cryst. Growth 2008, 310, 4330–4335. [Google Scholar] [CrossRef]
  4. Lv, G.; Chen, D.; Yang, X.; Ma, W.; Luo, T.; Wei, K.; Zhou, Y.; Zheng, G. Numerical simulation and experimental verification of vacuum directional solidification process for multicrystalline silicon. Vacuum 2015, 116, 96–103. [Google Scholar] [CrossRef]
  5. Wang, K. Fundations of Materials Engineering, 2nd ed.; Tsinghua University Press: Beijing, China, 2009. [Google Scholar]
  6. Nguyen, T.H.; Liao, S.H.; Chen, J.C.; Chen, C.H.; Huang, Y.H.; Yang, C.J.; Lin, H.W.; Nguyen, H.B. Effects of the hot zone design during the growth of large size multi-crystalline silicon ingots by the seeded directional solidification process. J. Cryst. Growth 2016, 452, 27–34. [Google Scholar] [CrossRef]
  7. Chen, L.G.; Dai, B. Optimization of power consumption on silicon directional solidification system by using numerical simulations. J. Cryst. Growth 2012, 354, 86–92. [Google Scholar] [CrossRef]
  8. Stelian, C. Numerical modeling of carbon distribution and precipitation during directional solidification of photovoltaic silicon. Int. J. Heat Mass Transf. 2019, 145, 118775. [Google Scholar] [CrossRef]
  9. Cablea, M.; Zaidat, K.; Gagnoud, A.; Nouri, A.; Chichignoud, G.; Delannoy, Y. Multi-crystalline silicon solidification under controlled forced convection. J. Cryst. Growth 2015, 417, 44–50. [Google Scholar] [CrossRef]
  10. Rao, S.; Chen, X.H.; Zhang, F.; He, L.; Luo, Y.; Xiong, H.; Hu, Y.; Wang, F.; Song, B. Influence of modified bottom insulation on the seeded directional solidification process for high-performance multi-crystalline silicon. Vacuum 2020, 172, 108969. [Google Scholar] [CrossRef]
  11. Pokland, E.S.A.; Heinze, V.; Meier, D.; Pätzold, O.; Stelter, M. Growth of multicrystalline silicon in a cone-shaped crucible. J. Cryst. Growth 2015, 416, 1–7. [Google Scholar]
  12. Lin, T.K.; Lin, C.H.; Chen, C.Y. Numerical analysis of steady and transient processes in a directional solidification system. Multi-Scale Multi-Phys. Mech. 2016, 1, 113–125. [Google Scholar] [CrossRef]
  13. Wu, Z.; Zhong, G.; Zhou, X.; Zhang, Z.; Wang, Z.; Chen, W.; Huang, X. Upgrade of the hot zone for large-size high performance multi-crystalline silicon ingot casting. J. Cryst. Growth 2016, 441, 58–63. [Google Scholar] [CrossRef]
  14. Srinivasan, M.; Nagarajan, S.G.; Ramasamy, P. Computational Study of Heat Transfer on Molten Silicon during Directional Solidification for Solar Cell Applications. Procedia Eng. 2015, 127, 1250–1255. [Google Scholar] [CrossRef] [Green Version]
  15. Ma, X.; Zheng, L.; Zhang, H.; Zhao, B.; Wang, C.; Xu, F. Thermal system design and optimization of an industrial silicon directional solidification system. J. Cryst. Growth 2011, 318, 288–292. [Google Scholar] [CrossRef]
  16. Trempa, C.R.M.; Jung, T.; Friedrich, J.; Müller, G. Modeling of incorporation of O, N, C and formation of related precipitates during directional solidification of silicon under consideration of variable processing parameters. J. Cryst. Growth 2010, 312, 878–885. [Google Scholar] [CrossRef]
  17. Yu, Q.; Liu, L.; Li, Z.; Shao, Y. Parameter study of traveling magnetic field for control of melt convection in directional solidification of crystalline silicon ingots. Int. J. Heat Fluid Flow 2018, 71, 55–67. [Google Scholar] [CrossRef]
  18. Li, Z.; Qi, X.; Liu, L.; Zhou, G. Numerical study of melt flow under the influence of heater-generating magnetic field during directional solidification of silicon ingots. J. Cryst. Growth 2018, 484, 78–85. [Google Scholar] [CrossRef]
  19. Kesavan, V.; Srinivasan, M.; Ramasamy, P. Numerical investigation of Directional Solidification process for improving multi-crystalline silicon ingot quality for photovoltaic applications. Mater. Lett. 2019, 241, 180–183. [Google Scholar] [CrossRef]
  20. Wu, B.; Stoddard, N.; Ma, R.H.; Clark, R. Bulk multicrystalline silicon growth for photovoltaic (PV) application. J. Cryst. Growth 2008, 310, 2178–2184. [Google Scholar] [CrossRef]
  21. Xie, G.; Lv, G.; Wang, Y.; Ma, W.; Yang, X.; Lei, Y. The influence of Marangoni effect on the growth quality of multi-crystalline silicon during the vacuum directional solidification process. Mater. Sci. Semicond. Process. 2019, 91, 124–132. [Google Scholar] [CrossRef]
  22. Karuppasamya, P.; Srinivasana, M.; Aravintha, K.; Ramasamya, P. Numerical Modelling on Modified Directional Solidification Process of Multi-crystalline Silicon Growth for Photovoltaic Applications. Mater. Today Proc. 2018, 5, 23014–23021. [Google Scholar] [CrossRef]
Figure 1. Two-dimensional representation of an axisymmetric directional solidification furnace for multicrystal Si and computational grids.
Figure 1. Two-dimensional representation of an axisymmetric directional solidification furnace for multicrystal Si and computational grids.
Crystals 10 01053 g001
Figure 2. Measurement and simulation temperature. (a) In the furnace without crucible; (b) in the operation condition of the ingot growth process.
Figure 2. Measurement and simulation temperature. (a) In the furnace without crucible; (b) in the operation condition of the ingot growth process.
Crystals 10 01053 g002
Figure 3. Temperature distribution and thermal stresses for the different power ratios at the half stage: (a) 1:1, (b) 1:2, (c) 1:3, and (d) 1:4.
Figure 3. Temperature distribution and thermal stresses for the different power ratios at the half stage: (a) 1:1, (b) 1:2, (c) 1:3, and (d) 1:4.
Crystals 10 01053 g003
Figure 4. Heat flux distributions in solids for different power ratios: (a) 1:1, (b) 1:2, (c) 1:3, and (d) 1:4.
Figure 4. Heat flux distributions in solids for different power ratios: (a) 1:1, (b) 1:2, (c) 1:3, and (d) 1:4.
Crystals 10 01053 g004
Figure 5. Radial and axial temperature distribution when the solidification fraction was 50%. (a) Radial temperature distribution on the horizontal plane (z = 1.0 m) and (b) axial temperature gradient between the plane of z = 1.0 m and the plane of z = 1.2 m.
Figure 5. Radial and axial temperature distribution when the solidification fraction was 50%. (a) Radial temperature distribution on the horizontal plane (z = 1.0 m) and (b) axial temperature gradient between the plane of z = 1.0 m and the plane of z = 1.2 m.
Crystals 10 01053 g005
Figure 6. (ad) Temperature distribution and thermal stress of different cases at the end state.
Figure 6. (ad) Temperature distribution and thermal stress of different cases at the end state.
Crystals 10 01053 g006
Figure 7. Diagram of the radiation model.
Figure 7. Diagram of the radiation model.
Crystals 10 01053 g007
Figure 8. Radiation heat flux and temperature distribution on the bottom surface.
Figure 8. Radiation heat flux and temperature distribution on the bottom surface.
Crystals 10 01053 g008
Figure 9. Positions of the block and heat flux for the furnace using different block positions: e = horizontal block, f = a 45°-angle-inclined block, and g = a vertical block.
Figure 9. Positions of the block and heat flux for the furnace using different block positions: e = horizontal block, f = a 45°-angle-inclined block, and g = a vertical block.
Crystals 10 01053 g009
Figure 10. (ad) Temperature and thermal stress distribution for different block positions at the half stage.
Figure 10. (ad) Temperature and thermal stress distribution for different block positions at the half stage.
Crystals 10 01053 g010
Figure 11. Radial and axial temperature distribution when the solidification fraction is 50%. (a) Radial temperature distribution on the horizontal plane (z = 1.0 m) and (b) axial temperature gradient between the plane of z = 1.0 m and the plane of z = 1.2 m.
Figure 11. Radial and axial temperature distribution when the solidification fraction is 50%. (a) Radial temperature distribution on the horizontal plane (z = 1.0 m) and (b) axial temperature gradient between the plane of z = 1.0 m and the plane of z = 1.2 m.
Crystals 10 01053 g011
Figure 12. (ad) Temperature distribution and thermal stress for different block positions at the end stage.
Figure 12. (ad) Temperature distribution and thermal stress for different block positions at the end stage.
Crystals 10 01053 g012
Table 1. Material properties of the furnace.
Table 1. Material properties of the furnace.
SiliconValue Value
Liquid specific heat capacity (J/kg∙K)710Liquid thermal conductivity (W/m∙K)22
Solid specific heat capacity (J/kg∙K)1000Solid thermal conductivity (W/m∙K)64
Emissivity (W/m∙K)0.3Latent heat of fusion (J/kg)1.587 × 106
Density of solid (kg/m3)2330Stress coefficient c11/c22/c331.653 × 1011
Density of liquid in 1600 K (kg/m3)2520Stress coefficient c12/c13/c236.393 × 1010
Phase change temperature (K)1687Stress coefficient c447.962 × 1010
CrucibleValue Value
Specific heat capacity (J/kg∙K)740Thermal conductivity (W/m∙K)4.8
Emissivity (W/m∙K)0.8
Heater/graphite/cooler plateValue Value
Specific heat capacity (J/kg∙K)740Thermal conductivity (W/m∙K)80
Emissivity (W/m∙K)0.8
InsulatorValue
Specific heat capacity (J/kg∙K)846Thermal conductivity (W/m∙K)0.4
Emissivity (W/m∙K)0.8
Table 2. Simulation cases.
Table 2. Simulation cases.
No.Heating Power and Ratio *Insulation Rising Velocity (m/s)
Case a26:26 kW, (1:1)2 × 10−6
Case b17.3:34.7 kW, (1:2)
Case c13:39 kW, (1:3)
Case d10.4:41.6 kW, (1:4)
* Note: the ratio presents the top power to side power.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Chen, C.; Liu, G.; Zhang, L.; Wang, G.; Hou, Y.; Li, Y. Numerical Simulation Study on the Front Shape and Thermal Stresses in Growing Multicrystalline Silicon Ingot: Process and Structural Design. Crystals 2020, 10, 1053. https://doi.org/10.3390/cryst10111053

AMA Style

Chen C, Liu G, Zhang L, Wang G, Hou Y, Li Y. Numerical Simulation Study on the Front Shape and Thermal Stresses in Growing Multicrystalline Silicon Ingot: Process and Structural Design. Crystals. 2020; 10(11):1053. https://doi.org/10.3390/cryst10111053

Chicago/Turabian Style

Chen, Chengmin, Guangxia Liu, Lei Zhang, Guodong Wang, Yanjin Hou, and Yan Li. 2020. "Numerical Simulation Study on the Front Shape and Thermal Stresses in Growing Multicrystalline Silicon Ingot: Process and Structural Design" Crystals 10, no. 11: 1053. https://doi.org/10.3390/cryst10111053

APA Style

Chen, C., Liu, G., Zhang, L., Wang, G., Hou, Y., & Li, Y. (2020). Numerical Simulation Study on the Front Shape and Thermal Stresses in Growing Multicrystalline Silicon Ingot: Process and Structural Design. Crystals, 10(11), 1053. https://doi.org/10.3390/cryst10111053

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop