- Email: [email protected]
Gas flow optimization during the cooling of multicrystalline silicon ingot
International Journal of Heat and Mass Transfer 84 (2015) 370–375
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Gas flow optimization during the cooling of multicrystalline silicon ingot S. Wang a, H.S. Fang a,⇑, C.J. Zhao a, Z. Zhang a, M.J. Zhang a, J.F. Xu b a b
School of Energy and Power Engineering, Huazhong University of Science and Technology, Wuhan 430074, China School of Mechanical Science and Engineering, Huazhong University of Science and Technology, Wuhan 430074, China
a r t i c l e
i n f o
Article history: Received 12 September 2013 Received in revised form 20 March 2014 Accepted 5 January 2015 Available online 22 January 2015 Keywords: Multicrystalline silicon Directional solidification Stresses Gas flow Cooling process
a b s t r a c t Multicrystalline silicon (mc-Si) produced by unidirectional solidification system is a crucial photovoltaic material due to its relatively high conversion efficiency and low cost. Defects related to thermal stresses, for example, dislocation, significantly affect the performance of the material. In the paper, a global transient model is applied to examine effects of gas flow on stress levels inner the silicon ingot during the cooling process. The maximum von Mises stresses under different inlet gas velocities are presented as a function of the cooling time. Stress level with a high inlet velocity at the initial cooling is slightly lower, but is much larger at the late cooling than that with a slow velocity. An optimized condition with variable gas velocities at the inlet is proposed to improve the quality of silicon ingot by reducing the stress level during the cooling process. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Recent years have witnessed impressive progress of the production of multicrystalline silicon (mc-Si) by the unidirectional solidification technique. However, silicon grown by this method has a relatively high dislocation density that reduces the minority carrier lifetime [1] and lowers photovoltaic conversion efficiency. Thermal stresses caused by inhomogeneous temperature field play an important role in generation and multiplication of dislocation during the solidification and cooling processes [2,3]. According to the theory [4], dislocations are easily generated inner the ingot where the effective stress exceeds the critical resolved shear stress (CRSS). Many researches on the stresses in silicon ingot have been carried out [5–13]. Chen et al. [5] used a transient and global model to study the thermal stress distribution in the silicon ingot for different solidification time, as well as the effect of crucible constrains. They suggested that a longer solidification time and lower constrains should be adopted for growing a silicon ingot with low thermal stress and low dislocation density. Takahashi et al. [6] investigated the generation mechanism of dislocation, and found that dislocations come into being at grain boundaries and propagated as crystal growth proceeds. Moreover, the shear stress on the slip plane around the grain boundary is likely the cause of dislocation formation. Fang et al. [7] performed parametric studies to ⇑ Corresponding author. Tel.: +86 27 87542618. E-mail addresses: [email protected] (H.S. Fang), [email protected] (J.F. Xu). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.01.035 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.
discuss the effect of furnace design on the interface shape and on the maximum von Mises stress. They concluded that the decrease of side insulation thickness, gas flow rate and distance from the bottom insulation to the heat exchanger block, as well as the increase of the top insulation thickness help to achieve a flat or slightly convex interface and the reduction of the maximum stress level and dislocation density in the growing ingot. Chen et al. [8] further studied the effect of thermal conductivity of a crucible on thermal stress and dislocations during the solidification using a three-dimensional global model, and concluded that temperature gradient in a silicon ingot as well as the melt-crystal interface shape should be controlled to reduce thermal stress and dislocations in a silicon ingot. Nakano et al. [9] considered the influence of the cooling rate on the dislocation density in mc-Si ingot during the unidirectional solidification process. They indicated that the maximum value of dislocation density decreases and that of residual stress increases with a fast cooling process. Zhou et al. [10] proposed a simplified model of dislocation density, and improved the traditional cooling process in the directional solidification of mc-Si. Their results showed that the final dislocation density could be reduced to about one fifth of that formed in normal cooling process. Wang et al. [11] built up a global transient model to study the effect of bottom insulation on thermal stress distribution in the mc-Si ingot during the cooling process, and found that the maximum von Mises stress is affected significantly by the insulation motion. In addition, generation and multiplication mainly occurs at the early cooling stage, and the lower part of the ingot always has no excessive stress, which means the region may have the best crystal quality. Zhao et al. [13] employed three methods to
S. Wang et al. / International Journal of Heat and Mass Transfer 84 (2015) 370–375
perform a quality assessment of silicon ingots produced in a directional solidification furnace. They concluded that the evaluation method based on thermo-elastic theory is preferred for efficient calculation in solidification process, and that for the whole solidification and annealing processes, the thermo-plastic theory and the thermo-creep theory are more suitable. They also found that large dislocation density is generated in the ingot during the annealing process. Although investigations on thermal stress distribution in multicrystalline silicon ingot have been conducted, there are few reports regarding the effect of gas velocity on the thermal stress during the cooling process. The main objective of the paper is to analyze the effect of the different gas velocities at the inlet on the cooling process, to present the relationship between the maximum von Mises stress and the gas velocities, and to propose an optimized gas inlet profile for the improvement of the crystal quality by reducing thermal-induced stress.
371
Fig. 2. The applied temperature profile on the heater and the motion of the bottom insulation during the cooling process.
2. Problem description and mathematical models 2.1. System configuration Schematic diagram of a typical directional solidification system is presented in Fig. 1. During the cooling process, the system components remain stationary except the bottom insulation that moves down directionally with a corresponding ramping-down of the heater power. To simplify the calculation, the heater power change is represented by a nonlinear temperature profile as described in Ref. [11]. The dimension of silicon ingot is 275 mm 880 mm 880 mm. The furnace information can also be found in Refs. [7,12]. The applied temperature profile along the heater and the motion of the bottom insulation are presented in Fig. 2. Thermal stresses generate during the solidification, release partially through the annealing, and pile up again during the cooling. The minimization of stress-related defects by reducing thermal stress has been an effective way to improve the crystal quality. In the actually industrial production, the silicon ingot is cuboids, and a three-dimensional model is required for an accurate analysis. However, the calculation time exceeds acceptance for an efficient study with the three-dimensionally transient model coupled with thermal field, flow field, stress field, and the tracking of the moving boundary due to the motion of the bottom insulation. Hence, a two-dimensional model, as applied in the Refs. [7,11], is adopted to save the computational load for a transient analysis of the cooling process.
2.2. Mathematical models In order to get the distribution of thermal stress inner the silicon ingot, energy conservation equation is solved by coupling with the incompressible N–S equation. Gas velocity at the inlet is less than 0.6 m/s, which makes a maximum value of Reynolds number less than 440. Therefore, the gas flow can be treated as laminar flow. The following assumptions are further applied: (1) all the components in the chamber are diffuse-gray body; (2) the Boussinesq approximation is used to solve natural convection. In order to make the approximation applicable, the inlet can be set at the top insulation, where the gas is heated up, to limit the temperature difference of the gas; (3) the gravity is negligible compared to the thermal stress in the stress balance equations. The governing equations can be written as
r~ u ¼ 0;
ð1Þ
@~ u q þ qð~ u rÞ~ u ¼ qrðmr ~ uÞ rp q0~ g bT ðT T 0 Þ; @t
ð2Þ
u is velocity vector, and q0 is the reference density; where t is time, ~ ~ g is the gravitational acceleration vector, bT is thermal expansion coefficient of the fluid, and T 0 is the reference temperature q and m are density and kinetic viscosity of the fluid, respectively. Temperature field is achieved by solving the energy conservative equation as
qi C p;i
* @T þ qi C p;i r ðT uÞ ¼ qi C p;i r ðai rTÞ þ r q00r;i ; @t
ð3Þ
where the subscript i indicates the components of the furnace; ai and Cp; i represent thermal diffusivity and specific heat of the material at the ith domain, respectively; q00r;i is the radiation source term. Thermal stress is calculated by
ruu 1 @ @ ðrrrr Þ þ ðrrz Þ ¼ 0; r @r @z r 1 @ @ ðrrrz Þ þ ðrzz Þ ¼ 0; r @r @z
ð4Þ
where rrr , ruu , and rzz are normal stresses in the radial, azimuthal and axial directions, respectively, and rrz is shear stress. The stress– strain relationship is required to solve the above equation. For anisotropic thermo-elastic solid body it is taken as:
rij ¼ C ijkl ekl bij ðT T ref Þdkl ; Fig. 1. Schematic diagram of a directional solidification furnace.
ð5Þ
where rij , ekl and dkl are stress tensor, strain tensor, and Kronecker delta tensor, respectively. bij is thermal expansion coefficient, and
372
S. Wang et al. / International Journal of Heat and Mass Transfer 84 (2015) 370–375
C ijkl is the elastic constant tensor. The above formula can be used for multicrystalline silicon with large grain size as discussed in our previous study [7]. For two-dimensional geometry, the above stress– strain equation is given as [5]:
0
1
10
0
1
rrr C 11 C 12 C 13 0 err b11 ðT T ref Þ C B CB C B B ruu C B C 12 C 22 C 23 0 CB eu/ b22 ðT T ref Þ C C B CB C B ¼ C B CB C; B B rzz C B C 13 C 23 C 33 0 CB ezz b ðT T ref Þ C 33 A @ A@ A @ rrz erz 0 0 0 C 44
ð6Þ
@v ; @z
ð7Þ
err ¼
@u ; @r
u r
euu ¼ ; ezz ¼
erz ¼
@ v @u þ @r @z
where Cij is the elastic constant; u and v are displacement components in the radial and axial directions, respectively. The parameters of the equations are discussed in the Ref. [7]. The von Mises stress is calculated from the stress components. The boundary conditions are described in the following. Temperatures at the exterior surface of the chamber and at the gas inlet are set to room temperature, which is reasonable considering the efficient cooling by water. Input power of the heater is reasonably substituted by a specific temperature profile to improve the convergence of the calculation. The applied temperature profile as shown in Fig. 2 is obtained by data fitting of the measured data under normal operating procedures. In our former study, we compared the measured temperature and the predicted temperature at TC2, and found that they fit quite well during the cooling stage [11]. It is noted that natural convection of the gas before entering the chamber, i.e., the part from inlet 1 to the top insulation 2 as shown in Fig. 1, is not calculated. The condition of convective flux is added to the outlet. The inlet velocities normal to the boundary are 0.6 m/s, 0.4 m/s, 0.3 m/s, 0.22 m/s, 0.05 m/s for Case 1, Case 2, Case 3, Case 4 and Case 5, respectively. All interfaces between gas and solid are non-slip except on the moving insulation, which is treated as moving wall with a vertical speed of 9e06 m/s. No * * traction boundary (r n ¼ 0) is employed on the top and side boundaries of the ingot. The weight of the ingot is applied at the bottom of the ingot. At the axis, an axisymmetric boundary condition is applied for all variables. As the bottom-insulation moves down, grids around the moving insulation experience different stretches or compressions. A large number of distorted meshes may appear to degrade the precision and convergence of numerical computation. Therefore, it is very important to keep good mesh equality in every discrete time step. In the current simulation, the Laplace smoothing is selected to constraint the meshing on the boundaries. The grid system can be adapted to eliminate low-quality meshes with the moving-mesh interface.
3. Results and discussion The influences of gas inlet velocities on transient cooling process are investigated. The cooling process takes about 390 min according to the industrial experience. The predicted distributions of thermal field, flow field and stress field at four different cooling stages are presented in Fig. 3. Temperature distribution in the furnace is represented by isothermal curves on the left of the axis, gas flow is marked by sequence numbers on the right of the axis, and von Mises stress is depicted by contour lines inner the silicon ingot. Radiative heat exchanges at all diffuse surfaces in the whole chamber play a dominant role under the environment of high temperature at the early cooling stage. Heat flux inner the solid components is transferred by heat conduction. The influence of convective heat transfer get more and more remarkable as the cooling proceeds as the system temperature decreases. For the sil-
icon ingot, temperature reduction is mainly attributed to heat extraction from the surfaces (AB, AD, and DC as shown in Fig. 3a (1)). Fig. 3a and b present simulation results with the inlet gas velocities of 0.05 m/s and 0.6 m/s, respectively. At the beginning of the cooling process (see Fig. 3a (1) and b (1)), both cases have similar temperature distributions. The highest furnace temperature is close to the heater, and the lowest ingot temperature locates at the bottom of the ingot. Heat flux comes from the heater, passes through the insulations and heat exchange block, and finally loses into the environment. Thermal distributions above surface AB are obviously different in Fig. 3a (1) and b (1), especially around the axis region. The isothermal lines in Fig. 3b are more parallel to the axis than those in Fig. 3a due to the high gas flow from the inlet along the axis. Gas velocity decreases gradually away from the axis. It can be expected that the major heat transfer regimes between the heat exchanger and the bottom insulation are thermal radiation and heat conduction due to a weak velocity there. In Fig. 3a (1) and b (1), the regions with lower thermal stresses are all located between the lower part of the ingot, and the maximum von Mises stresses both appear at the corners of the ingot. However, the maximum stress in Fig. 3a (1) reaches up to about 20 MPa, which is slightly higher than that in Fig. 3b (1), due to the slower gas cooling. After one hundred minutes cooling as shown in Fig. 3a (2) and b (2), the gap between the heat exchanger and the bottom insulation gets larger. The cold gas from the inlet enhances convective heat transfer on the surface AB. Temperature distribution inner the furnace changes dramatically due to the drop of the heater temperature. The high temperature regions have changed from the heater to the ingot, and the pattern of isotherms change from parallel lines to concentric ellipses in the ingot, which indicates that heat flux transfers in all directions from the center of silicon ingot to the environment. One can further find that the center of the high temperature zone in Fig. 3a (2) is located on the axis, while in Fig. 3b (2) it is away from the axis. Such patterns of temperature distributions are affected by the cooling gas, which leads to specific stress fields in the ingots. The maximum von Mises stress is still located at the corner of the ingot, however, regions with minimum stresses move down slightly. Conclusively, a higher gas velocity usually makes the ingot cool down more rapid, and results in larger stress level in the ingot. As the cooling process proceeds to 250 min, the maximum temperatures in the ingot drop to 1211 K and 1204 K for the cases in Fig. 3a (3) and b (3), respectively. Large stresses appear at the upper surface of the ingot due to locally large temperature gradients caused by the enhanced convective cooling. The maximum von Mises stresses move from the bottom corner of the ingot to the center of surface AB for both of the cases. Nevertheless, the variation of the maximum von Mises stress in Fig. 3a (3) and b (3) are different. For the slow gas velocity, the maximum value drops from about 22 MPa to 14 MPa in Fig. 3a (3) due to a slow reduction of temperature difference inner the silicon ingot. For the fast gas velocity, it goes up from about 21 MPa to 26 MPa in Fig. 3b (3), which may be caused by an increase of temperature gradient at surface AB due to a fast cooling rate. Fig. 3a (4) and b (4) illustrate temperature, velocity and von Mises stress distributions at the end of the cooling process. The bottom insulation reaches the lowest position, and the gap between the moving insulation block and the side insulation is the largest. Temperature at the centre of the silicon ingot in Fig. 3a (4) drops to 947 K, and that in Fig. 3b (4) is about 940 K. The pattern of thermal stress is similar to that at the cooling time of 250 min. Fig. 4 shows the maximum von Mises stress as a function of the inlet velocity. It can be found that at the beginning of the cooling stage the stress decreases sharply, which may be caused by a local
S. Wang et al. / International Journal of Heat and Mass Transfer 84 (2015) 370–375
(1) t = 25 min
373
(2) t = 100 min
(3) t = 250 min
(4) t = 393 min
(a)
(1) t = 25 min
(2) t = 100 min
(3) t = 250 min
(4) t = 393 min
(b) Fig. 3. Distributions of temperature (left), von Mises stress (in the silicon ingot), and velocity (right) during the cooling process at different time with a gas velocity of (a) 0.05 m/s, and (b) 0.6 m/s.
decrease of temperature gradients at the bottom of the ingot as illustrated in Fig. 3a (1) and b (1). During the cooling, the heater temperature goes down, and the temperature difference of the ingot and the heater, as a result, decrease significantly at the initial cooling stage. After a short cooling time, around 25 min later, the stress increases gradually. A peak value of about 22.5 MPa and
21.4 MPa is achieved for the gas velocity of 0.05 m/s and 0.6 m/s, respectively. Afterwards, the stress changes in a shape of ‘‘V’’. For the high velocities, the second peak can be obtained as marked by the green dots. For Case 1 with a slow gas velocity, the second peak does not present. The value at the second peak is the largest for the fastest gas flow (Case 1) with a value of about 26 MPa. The
374
S. Wang et al. / International Journal of Heat and Mass Transfer 84 (2015) 370–375
Fig. 4. The maximum von Mises stresses as a function of cooling time under constant velocities.
Fig. 6. The maximum von Mises stress as a function of cooling time for the different velocity profiles.
maximum von Mises stress of Case 1 is always larger than that of Case 5 before the cooling time of 150 min. Afterwards, the stress level of Case 5 is much larger than that of Case 1. It is because that the high gas velocity leads to rapid cooling. Therefore, the stress level builds up in the ingot gradually. To minimize thermalinduced stress during the cooling process definitely is necessary to reduce stress-related defects in the ingot. Conclusively, the maximum von Mises stress increases with gas velocity after 200 min cooling, which means that the slow cooling rate takes less time to reach the stress level as the fast cooling does. Specially, when the stress value reaches 16 MPa that is believed to be safe for silicon, it takes 200 min for Case 1, 275 min for Case 2, 375 min for Case 3, and more time for Case 4–5. The cooling stage can be stopped when the maximum stress is below the value. Hence, from the viewpoint of stress level, a slow cooling rate is expected to save cooling time. However, from the point of the final system temperature, it is apparent that slowing down gas rate will increase the total cooling time. From the above analysis, one may wonder what happens if the gas inlet is high at the initial cooling stage, and then becomes slow at the late cooling stage. There may be an optimum velocity minimizing the maximum von Mises stress at the whole cooling process.
To verify the thought, a complete study of various inlet velocities is presented. Four velocity profiles are discussed in Fig. 5. The starting and ending gas velocities are 0.6 m/s and 0.05 m/s, respectively, corresponding to Case 5 and Case 1 in Fig. 4. Gas velocity changes from 0.6 m/s to 0.05 m/s with different ramp-down rates as depicted by v1, v2, v3 and v4 in Fig. 5. Fig. 6 presents the evolution of the maximum von Mises stress as a function of the cooling time. The curves a, b, c and d, are related to the velocity profiles of v1, v2, v3 and v4, respectively. The shapes of the curves are similar to that of Case 1 in Fig. 4, when a slow gas velocity is given. The peak value ranges from 21.4 MPa to 22.2 MPa, and the minimum stress is about 13.7 MPa. It is noted that for most of the cooling time the stresses of the four profiles have close values, and that during the cooling time from 60 min to 200 min the differences of the stresses are obvious. In Fig. 7, the optimized curve almost overlaps upon that for a gas velocity of 0.6 m/s before the starting cooling of 150 min, and for a gas velocity of 0.05 m/s after the late cooling of 240 min. Hence, the velocity v3 with the lowest stress level can be taken as the optimal value. Although the differences of the stress values are small, the averaged stress level in the ingot could be improved greatly. For the regions with stress levels close to the critical value to generate defects, such an optimized profile is especially
Fig. 5. Gas velocities as a function of the cooling time.
Fig. 7. The maximum von Mises stress as a function of the cooling time.
S. Wang et al. / International Journal of Heat and Mass Transfer 84 (2015) 370–375
important. Fig. 7 redraws the maximum von Mises stresses of Case 1 (Fig. 4), Case 5 (Fig. 4) and curve c (Fig. 5). One can see clearly that the optimized condition of gas inlet is very effective to reduce thermal stress in the mc-Si ingot during the cooling process. Therefore, the quality of the silicon ingot can be well improved without an extension of the cooling time.
375
Acknowledgements The work is supported in part by the National Natural Science Foundation of China (No. 51476068), and in part by Key Foundation of Hubei Province (2013CFA080). References
4. Conclusions In the paper, a global and transient model is applied to study the cooling process of multicrystalline silicon in a directional solidification system. Temperature, velocity and thermal stress distributions at the specific cooling time are presented with different gas inlet velocities. The following conclusions can be drawn: (1) At the initial cooling stage, the region with the largest stress is at the bottom of the ingot. Then, it moves to the centre and side of the ingot. At the late cooling stage, the region close to the top surface of the ingot has the largest stress. (2) Thermal stress evolution in the ingot during the cooling is influenced by the inlet gas significantly. It decreases sharply at the initial cooling following by a gradual increase. Then, a ‘‘v’’ shape of the curve can be predicted for a fast gas velocity at the inlet. (3) A fast gas inlet causes a relatively lower stress at the early cooling and a very high stress at the late cooling; however, a slow gas inlet results in contrary results. Therefore, an optimized profile of the gas inlet velocity during the cooling process is proposed, during which a fast gas inlet at the initial cooling changes to a slow gas inlet at the late cooling by a certain ramping-down rate. The optimized profile is validated numerically to reduce the stress level during the completely cooling process without change of the cooling time. Conflict of interest None declared.
[1] K. Arafune, T. Sasaki, F. Wakabayashi, Y. Terada, Y. Ohshita, M. Yamaguchi, Study on defects and impurities in cast-grown polycrystalline silicon substrates for solar cells, Phys. B Condens. Matter 376–377 (2006) 236–239. [2] N. Miyazaki, Y. Kuroda, M. Sakaguchi, Dislocation density analyses of GaAs bulk single crystal during growth process (effects of crystal anisotropy), J. Cryst. Growth 218 (2000) 221–231. [3] X. Chen, S. Nakano, K. Kakimoto, Three-dimensional global analysis of thermal stress and dislocations in a silicon ingot during a unidirectional solidification process with a square crucible, J. Cryst. Growth 312 (2010) 3261–3266. [4] A.S. Jordan, R. Caruso, A. Von Neida, A thermoelastic analysis of dislocation generation in pulled GaAs crystals, Bell Syst. Tech. J. 59 (1980) 593–637. [5] X.J. Chen, S. Nakano, L.J. Liu, K. Kakimoto, Study on thermal stress in a silicon ingot during a unidirectional solidification process, J. Cryst. Growth 310 (2008) 4330–4335. [6] I. Takahashi, N. Usami, K. Kutsukake, G. Stokkan, K. Morishita, K. Nakajima, Generation mechanism of dislocations during directional solidification of multicrystalline silicon using artificially designed seed, J. Cryst. Growth 312 (2010) 897–901. [7] H.S. Fang, S. Wang, L. Zhou, N.G. Zhou, M.H. Lin, Influence of furnace design on the thermal stress during directional solidification of multicrystalline silicon, J. Cryst. Growth 346 (2012) 5–11. [8] X.J. Chen, S. Nakano, K. Kakimoto, 3D numerical analysis of the influence of material property of a crucible on stress and dislocation in multicrystalline silicon for solar cells, J. Cryst. Growth 318 (2011) 259–264. [9] S. Nakano, X.J. Chen, B. Gao, K. Kakimoto, Numerical analysis of cooling rate dependence on dislocation density in multicrystalline silicon for solar cells, J. Cryst. Growth 318 (2011) 280–282. [10] N. Zhou, M. Lin, L. Zhou, Q. Hu, H. Fang, S. Wang, A modified cooling process in directional solidification of multicrystalline silicon, J. Cryst. Growth 381 (2013) 22–26. [11] S. Wang, H.S. Fang, G.H. Wei, L. Zhou, N.G. Zhou, Transient analysis of the cooling process and influence of bottom insulation on the stress in the multicrystalline silicon ingot, Cryst. Res. Technol. 48 (2013) 454–463. [12] X. Ma, L. Zheng, H. Zhang, B. Zhao, C. Wang, F. Xu, Thermal system design and optimization of an industrial silicon directional solidification system, J. Cryst. Growth 318 (2011) 288–292. [13] Wenhan Zhao, Lijun Liu, Lei Sun, A’nan Geng, Quality evaluation of multicrystalline silicon ingots produced in a directional solidification furnace with different theories, J. Cryst. Growth 401 (2014) 296–301.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Gas flow optimization during the cooling of multicrystalline silicon ingot S. Wang a, H.S. Fang a,⇑, C.J. Zhao a, Z. Zhang a, M.J. Zhang a, J.F. Xu b a b
School of Energy and Power Engineering, Huazhong University of Science and Technology, Wuhan 430074, China School of Mechanical Science and Engineering, Huazhong University of Science and Technology, Wuhan 430074, China
a r t i c l e
i n f o
Article history: Received 12 September 2013 Received in revised form 20 March 2014 Accepted 5 January 2015 Available online 22 January 2015 Keywords: Multicrystalline silicon Directional solidification Stresses Gas flow Cooling process
a b s t r a c t Multicrystalline silicon (mc-Si) produced by unidirectional solidification system is a crucial photovoltaic material due to its relatively high conversion efficiency and low cost. Defects related to thermal stresses, for example, dislocation, significantly affect the performance of the material. In the paper, a global transient model is applied to examine effects of gas flow on stress levels inner the silicon ingot during the cooling process. The maximum von Mises stresses under different inlet gas velocities are presented as a function of the cooling time. Stress level with a high inlet velocity at the initial cooling is slightly lower, but is much larger at the late cooling than that with a slow velocity. An optimized condition with variable gas velocities at the inlet is proposed to improve the quality of silicon ingot by reducing the stress level during the cooling process. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Recent years have witnessed impressive progress of the production of multicrystalline silicon (mc-Si) by the unidirectional solidification technique. However, silicon grown by this method has a relatively high dislocation density that reduces the minority carrier lifetime [1] and lowers photovoltaic conversion efficiency. Thermal stresses caused by inhomogeneous temperature field play an important role in generation and multiplication of dislocation during the solidification and cooling processes [2,3]. According to the theory [4], dislocations are easily generated inner the ingot where the effective stress exceeds the critical resolved shear stress (CRSS). Many researches on the stresses in silicon ingot have been carried out [5–13]. Chen et al. [5] used a transient and global model to study the thermal stress distribution in the silicon ingot for different solidification time, as well as the effect of crucible constrains. They suggested that a longer solidification time and lower constrains should be adopted for growing a silicon ingot with low thermal stress and low dislocation density. Takahashi et al. [6] investigated the generation mechanism of dislocation, and found that dislocations come into being at grain boundaries and propagated as crystal growth proceeds. Moreover, the shear stress on the slip plane around the grain boundary is likely the cause of dislocation formation. Fang et al. [7] performed parametric studies to ⇑ Corresponding author. Tel.: +86 27 87542618. E-mail addresses: [email protected] (H.S. Fang), [email protected] (J.F. Xu). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.01.035 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.
discuss the effect of furnace design on the interface shape and on the maximum von Mises stress. They concluded that the decrease of side insulation thickness, gas flow rate and distance from the bottom insulation to the heat exchanger block, as well as the increase of the top insulation thickness help to achieve a flat or slightly convex interface and the reduction of the maximum stress level and dislocation density in the growing ingot. Chen et al. [8] further studied the effect of thermal conductivity of a crucible on thermal stress and dislocations during the solidification using a three-dimensional global model, and concluded that temperature gradient in a silicon ingot as well as the melt-crystal interface shape should be controlled to reduce thermal stress and dislocations in a silicon ingot. Nakano et al. [9] considered the influence of the cooling rate on the dislocation density in mc-Si ingot during the unidirectional solidification process. They indicated that the maximum value of dislocation density decreases and that of residual stress increases with a fast cooling process. Zhou et al. [10] proposed a simplified model of dislocation density, and improved the traditional cooling process in the directional solidification of mc-Si. Their results showed that the final dislocation density could be reduced to about one fifth of that formed in normal cooling process. Wang et al. [11] built up a global transient model to study the effect of bottom insulation on thermal stress distribution in the mc-Si ingot during the cooling process, and found that the maximum von Mises stress is affected significantly by the insulation motion. In addition, generation and multiplication mainly occurs at the early cooling stage, and the lower part of the ingot always has no excessive stress, which means the region may have the best crystal quality. Zhao et al. [13] employed three methods to
S. Wang et al. / International Journal of Heat and Mass Transfer 84 (2015) 370–375
perform a quality assessment of silicon ingots produced in a directional solidification furnace. They concluded that the evaluation method based on thermo-elastic theory is preferred for efficient calculation in solidification process, and that for the whole solidification and annealing processes, the thermo-plastic theory and the thermo-creep theory are more suitable. They also found that large dislocation density is generated in the ingot during the annealing process. Although investigations on thermal stress distribution in multicrystalline silicon ingot have been conducted, there are few reports regarding the effect of gas velocity on the thermal stress during the cooling process. The main objective of the paper is to analyze the effect of the different gas velocities at the inlet on the cooling process, to present the relationship between the maximum von Mises stress and the gas velocities, and to propose an optimized gas inlet profile for the improvement of the crystal quality by reducing thermal-induced stress.
371
Fig. 2. The applied temperature profile on the heater and the motion of the bottom insulation during the cooling process.
2. Problem description and mathematical models 2.1. System configuration Schematic diagram of a typical directional solidification system is presented in Fig. 1. During the cooling process, the system components remain stationary except the bottom insulation that moves down directionally with a corresponding ramping-down of the heater power. To simplify the calculation, the heater power change is represented by a nonlinear temperature profile as described in Ref. [11]. The dimension of silicon ingot is 275 mm 880 mm 880 mm. The furnace information can also be found in Refs. [7,12]. The applied temperature profile along the heater and the motion of the bottom insulation are presented in Fig. 2. Thermal stresses generate during the solidification, release partially through the annealing, and pile up again during the cooling. The minimization of stress-related defects by reducing thermal stress has been an effective way to improve the crystal quality. In the actually industrial production, the silicon ingot is cuboids, and a three-dimensional model is required for an accurate analysis. However, the calculation time exceeds acceptance for an efficient study with the three-dimensionally transient model coupled with thermal field, flow field, stress field, and the tracking of the moving boundary due to the motion of the bottom insulation. Hence, a two-dimensional model, as applied in the Refs. [7,11], is adopted to save the computational load for a transient analysis of the cooling process.
2.2. Mathematical models In order to get the distribution of thermal stress inner the silicon ingot, energy conservation equation is solved by coupling with the incompressible N–S equation. Gas velocity at the inlet is less than 0.6 m/s, which makes a maximum value of Reynolds number less than 440. Therefore, the gas flow can be treated as laminar flow. The following assumptions are further applied: (1) all the components in the chamber are diffuse-gray body; (2) the Boussinesq approximation is used to solve natural convection. In order to make the approximation applicable, the inlet can be set at the top insulation, where the gas is heated up, to limit the temperature difference of the gas; (3) the gravity is negligible compared to the thermal stress in the stress balance equations. The governing equations can be written as
r~ u ¼ 0;
ð1Þ
@~ u q þ qð~ u rÞ~ u ¼ qrðmr ~ uÞ rp q0~ g bT ðT T 0 Þ; @t
ð2Þ
u is velocity vector, and q0 is the reference density; where t is time, ~ ~ g is the gravitational acceleration vector, bT is thermal expansion coefficient of the fluid, and T 0 is the reference temperature q and m are density and kinetic viscosity of the fluid, respectively. Temperature field is achieved by solving the energy conservative equation as
qi C p;i
* @T þ qi C p;i r ðT uÞ ¼ qi C p;i r ðai rTÞ þ r q00r;i ; @t
ð3Þ
where the subscript i indicates the components of the furnace; ai and Cp; i represent thermal diffusivity and specific heat of the material at the ith domain, respectively; q00r;i is the radiation source term. Thermal stress is calculated by
ruu 1 @ @ ðrrrr Þ þ ðrrz Þ ¼ 0; r @r @z r 1 @ @ ðrrrz Þ þ ðrzz Þ ¼ 0; r @r @z
ð4Þ
where rrr , ruu , and rzz are normal stresses in the radial, azimuthal and axial directions, respectively, and rrz is shear stress. The stress– strain relationship is required to solve the above equation. For anisotropic thermo-elastic solid body it is taken as:
rij ¼ C ijkl ekl bij ðT T ref Þdkl ; Fig. 1. Schematic diagram of a directional solidification furnace.
ð5Þ
where rij , ekl and dkl are stress tensor, strain tensor, and Kronecker delta tensor, respectively. bij is thermal expansion coefficient, and
372
S. Wang et al. / International Journal of Heat and Mass Transfer 84 (2015) 370–375
C ijkl is the elastic constant tensor. The above formula can be used for multicrystalline silicon with large grain size as discussed in our previous study [7]. For two-dimensional geometry, the above stress– strain equation is given as [5]:
0
1
10
0
1
rrr C 11 C 12 C 13 0 err b11 ðT T ref Þ C B CB C B B ruu C B C 12 C 22 C 23 0 CB eu/ b22 ðT T ref Þ C C B CB C B ¼ C B CB C; B B rzz C B C 13 C 23 C 33 0 CB ezz b ðT T ref Þ C 33 A @ A@ A @ rrz erz 0 0 0 C 44
ð6Þ
@v ; @z
ð7Þ
err ¼
@u ; @r
u r
euu ¼ ; ezz ¼
erz ¼
@ v @u þ @r @z
where Cij is the elastic constant; u and v are displacement components in the radial and axial directions, respectively. The parameters of the equations are discussed in the Ref. [7]. The von Mises stress is calculated from the stress components. The boundary conditions are described in the following. Temperatures at the exterior surface of the chamber and at the gas inlet are set to room temperature, which is reasonable considering the efficient cooling by water. Input power of the heater is reasonably substituted by a specific temperature profile to improve the convergence of the calculation. The applied temperature profile as shown in Fig. 2 is obtained by data fitting of the measured data under normal operating procedures. In our former study, we compared the measured temperature and the predicted temperature at TC2, and found that they fit quite well during the cooling stage [11]. It is noted that natural convection of the gas before entering the chamber, i.e., the part from inlet 1 to the top insulation 2 as shown in Fig. 1, is not calculated. The condition of convective flux is added to the outlet. The inlet velocities normal to the boundary are 0.6 m/s, 0.4 m/s, 0.3 m/s, 0.22 m/s, 0.05 m/s for Case 1, Case 2, Case 3, Case 4 and Case 5, respectively. All interfaces between gas and solid are non-slip except on the moving insulation, which is treated as moving wall with a vertical speed of 9e06 m/s. No * * traction boundary (r n ¼ 0) is employed on the top and side boundaries of the ingot. The weight of the ingot is applied at the bottom of the ingot. At the axis, an axisymmetric boundary condition is applied for all variables. As the bottom-insulation moves down, grids around the moving insulation experience different stretches or compressions. A large number of distorted meshes may appear to degrade the precision and convergence of numerical computation. Therefore, it is very important to keep good mesh equality in every discrete time step. In the current simulation, the Laplace smoothing is selected to constraint the meshing on the boundaries. The grid system can be adapted to eliminate low-quality meshes with the moving-mesh interface.
3. Results and discussion The influences of gas inlet velocities on transient cooling process are investigated. The cooling process takes about 390 min according to the industrial experience. The predicted distributions of thermal field, flow field and stress field at four different cooling stages are presented in Fig. 3. Temperature distribution in the furnace is represented by isothermal curves on the left of the axis, gas flow is marked by sequence numbers on the right of the axis, and von Mises stress is depicted by contour lines inner the silicon ingot. Radiative heat exchanges at all diffuse surfaces in the whole chamber play a dominant role under the environment of high temperature at the early cooling stage. Heat flux inner the solid components is transferred by heat conduction. The influence of convective heat transfer get more and more remarkable as the cooling proceeds as the system temperature decreases. For the sil-
icon ingot, temperature reduction is mainly attributed to heat extraction from the surfaces (AB, AD, and DC as shown in Fig. 3a (1)). Fig. 3a and b present simulation results with the inlet gas velocities of 0.05 m/s and 0.6 m/s, respectively. At the beginning of the cooling process (see Fig. 3a (1) and b (1)), both cases have similar temperature distributions. The highest furnace temperature is close to the heater, and the lowest ingot temperature locates at the bottom of the ingot. Heat flux comes from the heater, passes through the insulations and heat exchange block, and finally loses into the environment. Thermal distributions above surface AB are obviously different in Fig. 3a (1) and b (1), especially around the axis region. The isothermal lines in Fig. 3b are more parallel to the axis than those in Fig. 3a due to the high gas flow from the inlet along the axis. Gas velocity decreases gradually away from the axis. It can be expected that the major heat transfer regimes between the heat exchanger and the bottom insulation are thermal radiation and heat conduction due to a weak velocity there. In Fig. 3a (1) and b (1), the regions with lower thermal stresses are all located between the lower part of the ingot, and the maximum von Mises stresses both appear at the corners of the ingot. However, the maximum stress in Fig. 3a (1) reaches up to about 20 MPa, which is slightly higher than that in Fig. 3b (1), due to the slower gas cooling. After one hundred minutes cooling as shown in Fig. 3a (2) and b (2), the gap between the heat exchanger and the bottom insulation gets larger. The cold gas from the inlet enhances convective heat transfer on the surface AB. Temperature distribution inner the furnace changes dramatically due to the drop of the heater temperature. The high temperature regions have changed from the heater to the ingot, and the pattern of isotherms change from parallel lines to concentric ellipses in the ingot, which indicates that heat flux transfers in all directions from the center of silicon ingot to the environment. One can further find that the center of the high temperature zone in Fig. 3a (2) is located on the axis, while in Fig. 3b (2) it is away from the axis. Such patterns of temperature distributions are affected by the cooling gas, which leads to specific stress fields in the ingots. The maximum von Mises stress is still located at the corner of the ingot, however, regions with minimum stresses move down slightly. Conclusively, a higher gas velocity usually makes the ingot cool down more rapid, and results in larger stress level in the ingot. As the cooling process proceeds to 250 min, the maximum temperatures in the ingot drop to 1211 K and 1204 K for the cases in Fig. 3a (3) and b (3), respectively. Large stresses appear at the upper surface of the ingot due to locally large temperature gradients caused by the enhanced convective cooling. The maximum von Mises stresses move from the bottom corner of the ingot to the center of surface AB for both of the cases. Nevertheless, the variation of the maximum von Mises stress in Fig. 3a (3) and b (3) are different. For the slow gas velocity, the maximum value drops from about 22 MPa to 14 MPa in Fig. 3a (3) due to a slow reduction of temperature difference inner the silicon ingot. For the fast gas velocity, it goes up from about 21 MPa to 26 MPa in Fig. 3b (3), which may be caused by an increase of temperature gradient at surface AB due to a fast cooling rate. Fig. 3a (4) and b (4) illustrate temperature, velocity and von Mises stress distributions at the end of the cooling process. The bottom insulation reaches the lowest position, and the gap between the moving insulation block and the side insulation is the largest. Temperature at the centre of the silicon ingot in Fig. 3a (4) drops to 947 K, and that in Fig. 3b (4) is about 940 K. The pattern of thermal stress is similar to that at the cooling time of 250 min. Fig. 4 shows the maximum von Mises stress as a function of the inlet velocity. It can be found that at the beginning of the cooling stage the stress decreases sharply, which may be caused by a local
S. Wang et al. / International Journal of Heat and Mass Transfer 84 (2015) 370–375
(1) t = 25 min
373
(2) t = 100 min
(3) t = 250 min
(4) t = 393 min
(a)
(1) t = 25 min
(2) t = 100 min
(3) t = 250 min
(4) t = 393 min
(b) Fig. 3. Distributions of temperature (left), von Mises stress (in the silicon ingot), and velocity (right) during the cooling process at different time with a gas velocity of (a) 0.05 m/s, and (b) 0.6 m/s.
decrease of temperature gradients at the bottom of the ingot as illustrated in Fig. 3a (1) and b (1). During the cooling, the heater temperature goes down, and the temperature difference of the ingot and the heater, as a result, decrease significantly at the initial cooling stage. After a short cooling time, around 25 min later, the stress increases gradually. A peak value of about 22.5 MPa and
21.4 MPa is achieved for the gas velocity of 0.05 m/s and 0.6 m/s, respectively. Afterwards, the stress changes in a shape of ‘‘V’’. For the high velocities, the second peak can be obtained as marked by the green dots. For Case 1 with a slow gas velocity, the second peak does not present. The value at the second peak is the largest for the fastest gas flow (Case 1) with a value of about 26 MPa. The
374
S. Wang et al. / International Journal of Heat and Mass Transfer 84 (2015) 370–375
Fig. 4. The maximum von Mises stresses as a function of cooling time under constant velocities.
Fig. 6. The maximum von Mises stress as a function of cooling time for the different velocity profiles.
maximum von Mises stress of Case 1 is always larger than that of Case 5 before the cooling time of 150 min. Afterwards, the stress level of Case 5 is much larger than that of Case 1. It is because that the high gas velocity leads to rapid cooling. Therefore, the stress level builds up in the ingot gradually. To minimize thermalinduced stress during the cooling process definitely is necessary to reduce stress-related defects in the ingot. Conclusively, the maximum von Mises stress increases with gas velocity after 200 min cooling, which means that the slow cooling rate takes less time to reach the stress level as the fast cooling does. Specially, when the stress value reaches 16 MPa that is believed to be safe for silicon, it takes 200 min for Case 1, 275 min for Case 2, 375 min for Case 3, and more time for Case 4–5. The cooling stage can be stopped when the maximum stress is below the value. Hence, from the viewpoint of stress level, a slow cooling rate is expected to save cooling time. However, from the point of the final system temperature, it is apparent that slowing down gas rate will increase the total cooling time. From the above analysis, one may wonder what happens if the gas inlet is high at the initial cooling stage, and then becomes slow at the late cooling stage. There may be an optimum velocity minimizing the maximum von Mises stress at the whole cooling process.
To verify the thought, a complete study of various inlet velocities is presented. Four velocity profiles are discussed in Fig. 5. The starting and ending gas velocities are 0.6 m/s and 0.05 m/s, respectively, corresponding to Case 5 and Case 1 in Fig. 4. Gas velocity changes from 0.6 m/s to 0.05 m/s with different ramp-down rates as depicted by v1, v2, v3 and v4 in Fig. 5. Fig. 6 presents the evolution of the maximum von Mises stress as a function of the cooling time. The curves a, b, c and d, are related to the velocity profiles of v1, v2, v3 and v4, respectively. The shapes of the curves are similar to that of Case 1 in Fig. 4, when a slow gas velocity is given. The peak value ranges from 21.4 MPa to 22.2 MPa, and the minimum stress is about 13.7 MPa. It is noted that for most of the cooling time the stresses of the four profiles have close values, and that during the cooling time from 60 min to 200 min the differences of the stresses are obvious. In Fig. 7, the optimized curve almost overlaps upon that for a gas velocity of 0.6 m/s before the starting cooling of 150 min, and for a gas velocity of 0.05 m/s after the late cooling of 240 min. Hence, the velocity v3 with the lowest stress level can be taken as the optimal value. Although the differences of the stress values are small, the averaged stress level in the ingot could be improved greatly. For the regions with stress levels close to the critical value to generate defects, such an optimized profile is especially
Fig. 5. Gas velocities as a function of the cooling time.
Fig. 7. The maximum von Mises stress as a function of the cooling time.
S. Wang et al. / International Journal of Heat and Mass Transfer 84 (2015) 370–375
important. Fig. 7 redraws the maximum von Mises stresses of Case 1 (Fig. 4), Case 5 (Fig. 4) and curve c (Fig. 5). One can see clearly that the optimized condition of gas inlet is very effective to reduce thermal stress in the mc-Si ingot during the cooling process. Therefore, the quality of the silicon ingot can be well improved without an extension of the cooling time.
375
Acknowledgements The work is supported in part by the National Natural Science Foundation of China (No. 51476068), and in part by Key Foundation of Hubei Province (2013CFA080). References
4. Conclusions In the paper, a global and transient model is applied to study the cooling process of multicrystalline silicon in a directional solidification system. Temperature, velocity and thermal stress distributions at the specific cooling time are presented with different gas inlet velocities. The following conclusions can be drawn: (1) At the initial cooling stage, the region with the largest stress is at the bottom of the ingot. Then, it moves to the centre and side of the ingot. At the late cooling stage, the region close to the top surface of the ingot has the largest stress. (2) Thermal stress evolution in the ingot during the cooling is influenced by the inlet gas significantly. It decreases sharply at the initial cooling following by a gradual increase. Then, a ‘‘v’’ shape of the curve can be predicted for a fast gas velocity at the inlet. (3) A fast gas inlet causes a relatively lower stress at the early cooling and a very high stress at the late cooling; however, a slow gas inlet results in contrary results. Therefore, an optimized profile of the gas inlet velocity during the cooling process is proposed, during which a fast gas inlet at the initial cooling changes to a slow gas inlet at the late cooling by a certain ramping-down rate. The optimized profile is validated numerically to reduce the stress level during the completely cooling process without change of the cooling time. Conflict of interest None declared.
[1] K. Arafune, T. Sasaki, F. Wakabayashi, Y. Terada, Y. Ohshita, M. Yamaguchi, Study on defects and impurities in cast-grown polycrystalline silicon substrates for solar cells, Phys. B Condens. Matter 376–377 (2006) 236–239. [2] N. Miyazaki, Y. Kuroda, M. Sakaguchi, Dislocation density analyses of GaAs bulk single crystal during growth process (effects of crystal anisotropy), J. Cryst. Growth 218 (2000) 221–231. [3] X. Chen, S. Nakano, K. Kakimoto, Three-dimensional global analysis of thermal stress and dislocations in a silicon ingot during a unidirectional solidification process with a square crucible, J. Cryst. Growth 312 (2010) 3261–3266. [4] A.S. Jordan, R. Caruso, A. Von Neida, A thermoelastic analysis of dislocation generation in pulled GaAs crystals, Bell Syst. Tech. J. 59 (1980) 593–637. [5] X.J. Chen, S. Nakano, L.J. Liu, K. Kakimoto, Study on thermal stress in a silicon ingot during a unidirectional solidification process, J. Cryst. Growth 310 (2008) 4330–4335. [6] I. Takahashi, N. Usami, K. Kutsukake, G. Stokkan, K. Morishita, K. Nakajima, Generation mechanism of dislocations during directional solidification of multicrystalline silicon using artificially designed seed, J. Cryst. Growth 312 (2010) 897–901. [7] H.S. Fang, S. Wang, L. Zhou, N.G. Zhou, M.H. Lin, Influence of furnace design on the thermal stress during directional solidification of multicrystalline silicon, J. Cryst. Growth 346 (2012) 5–11. [8] X.J. Chen, S. Nakano, K. Kakimoto, 3D numerical analysis of the influence of material property of a crucible on stress and dislocation in multicrystalline silicon for solar cells, J. Cryst. Growth 318 (2011) 259–264. [9] S. Nakano, X.J. Chen, B. Gao, K. Kakimoto, Numerical analysis of cooling rate dependence on dislocation density in multicrystalline silicon for solar cells, J. Cryst. Growth 318 (2011) 280–282. [10] N. Zhou, M. Lin, L. Zhou, Q. Hu, H. Fang, S. Wang, A modified cooling process in directional solidification of multicrystalline silicon, J. Cryst. Growth 381 (2013) 22–26. [11] S. Wang, H.S. Fang, G.H. Wei, L. Zhou, N.G. Zhou, Transient analysis of the cooling process and influence of bottom insulation on the stress in the multicrystalline silicon ingot, Cryst. Res. Technol. 48 (2013) 454–463. [12] X. Ma, L. Zheng, H. Zhang, B. Zhao, C. Wang, F. Xu, Thermal system design and optimization of an industrial silicon directional solidification system, J. Cryst. Growth 318 (2011) 288–292. [13] Wenhan Zhao, Lijun Liu, Lei Sun, A’nan Geng, Quality evaluation of multicrystalline silicon ingots produced in a directional solidification furnace with different theories, J. Cryst. Growth 401 (2014) 296–301.