Numerical Analysis of the Dislocation Density in Multicrystalline Silicon for Solar Cells by the Vertical Bridgman Process

Abstract

We studied the effects of cooling process on the generation of dislocations in multicrystalline silicon grown by the vertical Bridgman process. From the temperature field obtained by a global model, the stress relaxation and multiplication of dislocations were calculated using the Haasen-Alexander-Sumino model. It was found that the multiplication of dislocations is higher in fast cooling processes. It was confirmed that residual stress is low at high temperatures because the movement of the dislocations relaxes the thermal strain, while the residual stress increases with decreasing temperature, because of reduced motion of dislocations and formation of a strain field at lower temperatures.

1. Introduction

Photovoltaic (PV) market installations around the world reached a high record of 27.4 GW in 2011, representing growth of 140% over the previous year [1]. Multicrystalline silicon (mc-Si) is the most widely used material for solar cells. Mc-Si is generally produced using a directional solidification process because of the mass productivity and cost-effectiveness of the process. However, it is known that a high dislocation density is commonly generated in the silicon ingot during this production process, which reduces minority carrier lifetime [2]. In addition, the large residual stress in the crystal causes large fracture rates from the mechanical yield perspective and an accumulation of impurities that affects the electrical performance [3]. It is therefore necessary to understand the variations in the dislocation density and the residual stress and reduce them to improve the conversion efficiency of the solar cells.

The multiplication of dislocations is mainly caused by thermal stress. Thermal stresses are generated in the silicon ingot by nonuniform temperature distribution, which is mainly affected by furnace design [4] and cooling process design [5]. When the effective stress exceeds the critical resolved shear stress (CRSS) of the material, dislocations multiply easily in the crystal. Numerical simulations of dislocations in mc-Si grown by the conventional casting method have been reported in the literature [6–8]. However, there have been very few reports on silicon crystals grown by the vertical Bridgman process [9]. Therefore, we have performed a numerical simulation study of the variation of the dislocation density and the residual stress in mc-silicon grown by the vertical Bridgman process. In particular, we focused on the effects of the cooling rate on the dislocation density and the residual stress.

2. Simulation Model and Computation Method

A global model [10, 11] was used to study the temperature distributions across the entire furnace for the unidirectional solidification and cooling processes. An axisymmetric furnace is shown in Figure 1. The growth system consists of the silicon melt, the silicon crystal, the crucible, the pedestal, the heater shield, and the heaters. We used one top heater and two side heaters, represented by heaters 12 and 13 in the figure, respectively. The heater power ratio of the top heater to the side heater was 26 : 74. In the vertical Bridgman process, we moved the crucibles and pedestals downward to solidify the melt slowly from the bottom to the top, as shown in Figure 2. This is because the heat shield (component 9) in Figure 1 prevents heat radiation from the top and side heaters, and the bottom parts of the furnace are heated less than the top parts of the furnace to enhance crystallization.

Figure 1 

Configuration inside the furnace (1: melt; 2: crystal; 3-4: crucibles; 5-6: pedestals; 7–11: heat shields; 12–15: heaters).

(a)

(b)

Figure 2 

Position of the crucible (a) at the beginning of solidification and (b) at the end of solidification.

We considered conductive heat transfer in all the furnace components, heat radiation exchange between all the diffuse surfaces in the furnace, and heat convection for melt flow during heat transfer calculations for the entire furnace. We coupled these elements and solved the temperature field in them by using a transient global model. In addition, the melt-solid interface shape was obtained by a dynamic tracking method [12]. The major assumptions made in our model were the following: (1) the geometry of the furnace configuration is axisymmetric, (2) the heat radiation is modeled as gray-body radiation, (3) the effect of natural gas flow on heat transfer is neglected in the present furnace, and (4) the melt flow in the crucible is laminar and incompressible. The governing equations for flow, conductive heat transfer, and heat radiation in the furnace under these assumptions can thus be expressed as follows: ∇⋅𝑉^→𝜌𝜕=0,_→𝑉+𝜕𝑡^→𝑉⋅∇𝑉^→=−∇𝑝+∇⋅𝜇∇𝑉^→+∇𝑉^→𝑇−𝜌𝑔^→𝛽_𝑇𝑇−𝑇_𝑚,𝜌𝑐𝜕𝑇+𝜕𝑡^→𝑉⋅∇𝑇=∇⋅(𝑘∇𝑇)+𝑆_𝑄,𝑞𝑥^→𝜀𝑥^→−_{→𝑥∗∈𝜕𝑉}𝐾𝑥^→,𝑥^→∗1−𝜀𝑥^→∗𝜀𝑥^→∗𝑞𝑥^→∗𝑑𝑆^∗=𝜎_𝑠𝑇𝑥^4→−_{→𝑥∗∈𝜕𝑉}𝐾𝑥^→,𝑥^→∗𝜎_𝑠𝑇𝑥^4→𝑑𝑆^∗, (1)

where 𝑉^→, 𝜌, 𝑝, 𝜇, 𝑔^→, 𝛽_𝑇, 𝑐, and 𝑘 represent the melt velocity, the melt density, the melt pressure, the melt viscosity, gravitational acceleration, the thermal expansion coefficient, the heat capacity, and the thermal conductivity, respectively. 𝑞(𝑥^→), 𝜀(𝑥^→), and 𝜎_𝑠 are the heat flux, the emissivity, and the Stefan-Boltzmann constant, respectively. 𝑆_𝑄 is the heat capacity per unit volume in the crucible. 𝑥^→ and 𝑥^→∗ are infinitesimal radiation surface elements on 𝜕𝑉. 𝑑𝑆^∗ is the area of the infinitesimal surface element 𝑥^→. 𝑘(𝑥^→,𝑥^→∗) is the surface view factor between𝑥^→and𝑥^→∗, for which full details are given in [12].

Based on the temperature distributions in the crystal, stress relaxation, time-dependent creep deformation, and the multiplication of the dislocations were calculated from the beginning of the solidification process to the end of the cooling process at room temperature. We assumed that (1) the silicon crystal is isotropic, and the dislocation density analysis is carried out only for the [001] direction of the silicon crystal; (2) the boundary conditions at the top and bottom surfaces of crystal are treated as free boundaries, and the boundary conditions at the ingot-wall interfaces are also treated as free boundaries, that is, 𝜎^→⋅^→𝑛=0, because some coated materials on the crucible walls are used to reduce the crucible constraints. In the dislocation model, the total strain is given by 𝜀̇_𝑖𝑗=𝜀̇^𝑒_𝑖𝑗+𝜀̇^𝑇_𝑖𝑗+𝜀̇^𝑐_𝑖𝑗, (2)

where 𝜀^𝑒_𝑖𝑗, 𝜀^𝑇_𝑖𝑗, and 𝜀^𝑐_𝑖𝑗 are the elastic strain, the thermal strain, and the creep strain, respectively. The creep strain is related to the dislocation density, which can be described using the Haasen-Alexander-Sumino (HAS) model [13, 14]. In the HAS model for a multiaxial stress state [15], the creep rate 𝜀̇^𝑐_𝑖𝑗 and the mobile dislocation density multiplication rate 𝑁̇_𝑚 are given by 𝜀̇^𝑐_𝑖𝑗=12𝑏𝑘₀𝜏_eff^𝑝−𝑄exp𝑁𝑘𝑇^𝑚1𝐽₂⋅𝑆_𝑖𝑗,𝑁̇_𝑚=𝐾𝑘₀𝜏_eff^𝑝+𝜆exp−𝑄𝑁𝑘𝑇^𝑚,𝜏_eff=𝐽₂−𝐷𝑁_𝑚,𝐽₂=1𝑆_𝑖𝑗𝑆_𝑖𝑗,𝑆_𝑖𝑗=𝜎_𝑖𝑗−𝜎_𝑘𝑘3𝛿_𝑖𝑗, (3)(4)(5)(6)(7)

where 𝑏 is the magnitude of the Burgers vector, 𝑘 is the Boltzmann constant, 𝑇 is the absolute temperature in the silicon crystal, 𝑄 is the Peierls potential, 𝑘₀, 𝐾, 𝑝, and 𝜆 are materials constants, 𝐷 is the strain hardening factor, and 𝑆_𝑖𝑗 and 𝐽₂ are the deviatoric stress and the second invariant of the deviatoric stress, respectively. The value of 𝜏_eff is set to equal zero when 𝜏_eff≤0, which means that the creep strain rate 𝜀̇^𝑐_𝑖𝑗 and the multiplication rate of the mobile dislocation density 𝑁̇_𝑚 are set to zero. The governing partial differential equations for momentum balance in the axisymmetric case can be written as 1𝑟𝜕𝜕𝑟𝑟𝜎_𝑟𝑟+𝜕𝜎𝜕𝑧^𝑟𝑧−𝜎_𝜃𝜃𝑟1=0,𝑟𝜕𝜕𝑟𝑟𝜎_𝑟𝑧+𝜕𝜎𝜕𝑧^𝑧𝑧=0, (8)

where 𝜎_𝑟𝑟, 𝜎_𝑧𝑧, 𝜎_𝜃𝜃, and 𝜎_𝑟𝑧 are the normal stresses in the radial, axial, and azimuthal directions and the shear stress, respectively. If the elastic strain of (2) is applied to Hooke’s law for an axisymmetric plastic body model, we obtain the stress relaxation as follows: 𝜎_𝑟𝑟𝜎_𝜃𝜃𝜎_𝑧𝑧𝜎_𝑟𝑧=𝑐₁₁𝑐₁₂𝑐₁₃0𝑐₁₂𝑐₂₂𝑐₂₃0𝑐₁₃𝑐₂₃𝑐₃₃0000𝑐₄₄𝜀_𝑟𝑟−𝜀^𝑇_𝑟𝑟−𝜀^𝑐_𝑟𝑟𝜀_𝜃𝜃−𝜀^𝑇_𝜃𝜃−𝜀^𝑐_𝜃𝜃𝜀_𝑧𝑧−𝜀^𝑇_𝑧𝑧−𝜀^𝑐_𝑧𝑧𝜀_𝑟𝑧−𝜀^𝑇_𝑟𝑧−𝜀^𝑐_𝑟𝑧, (9a)

where 𝑐_𝑖𝑗 is the elastic coefficient of silicon. For cubic crystal lattices, there are only three independent coefficients: 𝑐₁₁, 𝑐₁₂, and 𝑐₄₄. The other coefficients can be expressed as 𝑐₁₁=𝑐₂₂=𝑐₃₃=165.77 GPa, 𝑐₁₂=𝑐₁₃=𝑐₂₃=63.9 GPA, and𝑐₄₄=79.62 GPa [16]. The total strain,𝜀_𝑟𝑟, 𝜀_𝜃𝜃, 𝜀_𝑧𝑧, and 𝜀_𝑟𝑧, can be calculated by 𝜀_𝑟𝑟=𝜕𝑢𝜕𝑟,𝜀_𝜃𝜃=𝑢𝑟,𝜀_𝑧𝑧=𝜕𝑣𝜕𝑧,𝜀_𝑟𝑧=𝜕𝑢+𝜕𝑧𝜕𝑣.𝜕𝑟 (9b)

The displacement boundary conditions are derived from the stress boundary conditions, that is, 𝜎^→⋅^→𝑛=0, at the top, bottom, and side surfaces. During numerical iteration process, the stress boundary conditions should be always satisfied. At the axial line, the radial displacement 𝑢 is set to zero. The creep strains are first solved from (3) and (4), and then substituted into (8) to obtain total strains. The basic steps are the following.(1)Give temperature field, and obtain thermal strain 𝜀^𝑇_𝑖𝑗.(2)Give initial dislocation density 𝑁_𝑚^(𝑘), creep strain 𝜀̇^𝑐(𝑘)_𝑖𝑗, and stresses 𝜎^(𝑘)_𝑖𝑗.(3)Update dislocation density 𝑁_𝑚^(𝑘+1), and 𝜀̇^𝑐(𝑘+1)_𝑖𝑗 by using (3) and (4) with 𝑁_𝑚^(𝑘), 𝜀̇^𝑐(𝑘)_𝑖𝑗, and 𝜎^(𝑘)_𝑖𝑗.(4)By substituting 𝜀̇^𝑐(𝑘+1)_𝑖𝑗 into (9a) and (9b) with given boundary conditions at all of the surfaces, obtain displacements 𝑢^(𝑘+1) and 𝑣^𝑘+1 and a new stress 𝜎^(𝑘+1)_𝑖𝑗.(5)If there is no convergence, return to step (3) and iterate again.

Figure 3 shows the calculation conditions, which are histories of the heater power and the crucible position. Initially, we performed a slow cooling process by moving the crucible downward at a rate of 0.25 mm/min until the solidification was complete. Then, we stopped the crucible movement and changed the heater power to perform calculations for different cooling rate conditions. Figure 4 shows the calculation results for the average temperature histories in a sample crystal. In the present paper, the average value of one variable is defined by volume-weighted average inside the global region. Five different cooling rates denoted by 1, 2, 3, 4, and 5 for cooling down to 900 K were set as −8.4, −6.4, −5.0, −3.8, and −3.0 K/min, respectively. Below 900 K, the dislocation density almost completely stops increasing [8].

Figure 3 

Analysis conditions for heater power (solid lines) and crucible positions (dashed lines).

Figure 4 

Average crystal temperature histories.

3. Results and Discussion

3.1. Variation of the Dislocation Density during the Cooling Process

Figures 5 and 6 show the calculation results for the histories and distributions of dislocation density at room temperature, respectively. It was found that the multiplication of dislocation density increased at higher cooling rates. Then, the increases in the dislocation density almost entirely stopped below 900 K for all cooling rates. The results showed that the dislocation density of cooling rate 1 is much larger than those of the other rates. This is because multiplication of the moving dislocations is caused by plastic deformation due to thermal stress relaxation. Under the analysis conditions of rate 1, the heat flow from the bottom of the crucible was high because of the high cooling rate. Thus, the thermal stress for cooling rate 1 became larger than those of the other cooling conditions at the beginning of the cooling process, as shown in Figure 7. The dislocation multiplication became higher than that of the other rates as well. This means that the moving dislocation density increased to relax the thermal stress when the thermal stress conditions changed dramatically. This indicates that maintaining the temperature distribution during the cooling process by slow cooling is important for reduction of the dislocation density.

Figure 5 

Average dislocation density versus average temperature during the cooling process.

(a)

(b)

(c)

Figure 6 

Dislocation density distributions at room temperature (1, 3, and 5).

Figure 7 

Average residual stress of plastic body (solid lines) and elastic body (dashed lines) versus the average temperature during the cooling process.

3.2. Variation of the Residual Stress during the Cooling Process

Figures 7 and 8 show the calculation results for the residual stress histories and the residual stress distributions at room temperature, respectively. To represent the residual stress level, the von Mises stress is introduced as 𝜎_von=32𝑆_𝑖𝑗𝑆_𝑖𝑗^1/2. (10)

(a)

(b)

(c)

Figure 8 

Residual stress distributions at room temperature (1, 3, and 5).

We calculated the average residual stresses based on the plastic and elastic body models, as shown in Figure 7. The elastic body model does not contain the effect of the moving dislocations; thus, during cooling process, when the temperature gradient gradually decreases inside crystal, the thermal stress or residual stress also gradually decreases. However, for viscoplastic body, the residual stress gradually increases even for decreased temperature gradient; the only explanation is that the increase of residual stress is due to the rapid increase of dislocation density (Figure 8).

In addition, the residual stress distributions at room temperature varied with the cooling rates. At the bottom of the crystal, higher cooling rates resulted in higher residual stress levels. In contrast, the residual stress of the elastic body model decreased with decreasing temperature and reached zero at room temperature at all cooling rates.

Figure 9 shows a schematic of results. In the elastic body model, the residual stress is high during cooling processes at high temperature because of the thermal strain. The residual stress then reaches zero at room temperature because the thermal strain is reduced to zero. However, the residual stress history of the plastic body is different from that of the elastic body. During cooling processes at high temperature, the thermal strain is relaxed because of the moving dislocations, and the residual stress is consequently low. We quantitatively confirmed that the amount of stress relaxation, which is equal to the difference in stress between the elastic and plastic bodies, increased because of the multiplication of dislocations as indicated by the yield point of the general stress-strain curve in Figure 10. Then, as the temperature decreases, the thermal activation effect decreases and the dislocation density motion remains unchanged. Consequently, the fixed dislocations produce a strain field and the residual stress is high at room temperature. Moreover, the results show that the residual stress for a high cooling rate at the bottom of a crystal became high because of the increased heat flow from the bottom of the crucible. In addition, the residual stress distributions at room temperature depend on the degree of reduction in thermal strain during the cooling process at high temperatures because of the multiplication of the dislocations. Thus, maintaining homogeneous temperature distributions in a crystal at high temperatures can reduce the residual stress at room temperature, because the dislocation density is uniformly distributed.

(a)

(b)

Figure 9 

Model of the relationship between the thermal strain and the residual stress. (a) Variation of the residual stress in the elastic body from high temperature to room temperature during cooling process. (b) Variation of the residual stress in the viscoplastic body from high temperature to room temperature during cooling process.

(a)

(b)

Figure 10 

The relationship between the relaxation stress and dislocations. (a) The model of the stress-strain curve at the yield point. (b) Comparison between the multiplication of the dislocation density and the relaxation stress in 1.

4. Conclusion

Based on the transient global model and the HAS model, we have studied the effects of the cooling process on the dislocation density in multicrystalline silicon grown by the vertical Bridgman process. We found that the dislocation multiplication becomes high at fast cooling rates because residual stress increases dramatically due to the multiplication of dislocations. Also, the residual stress is low at high temperatures because of the movement of dislocations. The residual stress then increases with decreasing temperature because the fixed dislocations produce a strain field. These results indicated that maintaining homogeneous temperature distributions reduces the dislocation density and the residual stress in the crystal.