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

Authors: Makoto Inoue; Satoshi Nakano; Hirofumi Harada; Yoshiji Miyamura; Bing Gao; Yoshihiro Kangawa; Koichi Kakimoto
Journal: International Journal of Photoenergy
 (2013)

---

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


---

## Body

## 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).Position of the crucible (a) at the beginning of solidification and (b) at the end of solidification.
(a)
(b)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:
(1)∇·V→=0,ρ(∂V→∂t+V→·∇V→)=-∇p+∇·[μ(∇V→+∇V→T)]-ρg→βT(T-Tm),ρc(∂T∂t+V→·∇T)=∇·(k∇T)+SQ,q(x→)ε(x→)-∫x→*∈∂VK(x→,x→*)1-ε(x→*)ε(x→*)q(x→*)dS*=σsT4(x→)-∫x→*∈∂VK(x→,x→*)σsT4(x→)dS*,whereV→, ρ, p, μ, g→, βT, c, and k 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. q(x→), ε(x→), and σs are the heat flux, the emissivity, and the Stefan-Boltzmann constant, respectively. SQ is the heat capacity per unit volume in the crucible. x→ and x→* are infinitesimal radiation surface elements on ∂V.  dS* is the area of the infinitesimal surface element  x→. k(x→,x→*) is the surface view factor betweenx→andx→*, 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, σ→·n→=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)ε˙ij=ε˙ije+ε˙ijT+ε˙ijc,whereεije, εijT, and εijc 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 ε˙ijc and the mobile dislocation density multiplication rate N˙m are given by
(3)ε˙ijc=12bk0(τeff)pexp(-QkT)Nm1J2·Sij,(4)N˙m=Kk0(τeff)p+λexp(-QkT)Nm,(5)τeff=J2-DNm,(6)J2=1SijSij,(7)Sij=σij-σkk3δij,whereb is the magnitude of the Burgers vector, k is the Boltzmann constant, T is the absolute temperature in the silicon crystal, Q is the Peierls potential, k0, K, p, and λ are materials constants, D is the strain hardening factor, and  Sij and J2 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 ε˙ijc and the multiplication rate of the mobile dislocation density N˙m are set to zero. The governing partial differential equations for momentum balance in the axisymmetric case can be written as
(8)1r∂∂r(rσrr)+∂∂z(σrz)-σθθr=0,1r∂∂r(rσrz)+∂∂z(σzz)=0,whereσrr, σzz, σθθ, and σrz 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:
(9a)(σrrσθθσzzσrz)=(c11c12c130c12c22c230c13c23c330000c44)(εrr-εrrT-εrrcεθθ-εθθT-εθθcεzz-εzzT-εzzcεrz-εrzT-εrzc),wherecij is the elastic coefficient of silicon. For cubic crystal lattices, there are only three independent coefficients: c11, c12, and c44. The other coefficients can be expressed as c11=c22=c33=165.77 GPa, c12=c13=c23=63.9 GPA, andc44=79.62 GPa  [16]. The total strain,εrr, εθθ, εzz, and εrz, can be calculated by
(9b)εrr=∂u∂r,εθθ=ur,εzz=∂v∂z,εrz=∂u∂z+∂v∂r.The displacement boundary conditions are derived from the stress boundary conditions, that is,σ→·n→=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 u 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εijT.(2)
Give initial dislocation densityNm(k), creep strain ε˙ijc(k), and stresses σij(k).(3)
Update dislocation densityNm(k+1), and ε˙ijc(k+1) by using (3) and (4) with Nm(k),  ε˙ijc(k), and σij(k).(4)
By substitutingε˙ijc(k+1) into (9a) and (9b) with given boundary conditions at all of the surfaces, obtain displacements  u(k+1) and vk+1 and a new stress σij(k+1).(5)
If  there is no convergence, return to step(3) and iterate again.Figure3 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

Figures5 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.Dislocation density distributions at room temperature (1, 3, and 5).
(a)
(b)
(c)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

Figures7 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
(10)σvon=(32SijSij)1/2.Residual stress distributions at room temperature (1, 3, and 5).
(a)
(b)
(c)We calculated the average residual stresses based on the plastic and elastic body models, as shown in Figure7. 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.Figure9 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.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)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.
(a)
(b)


## 3.1. Variation of the Dislocation Density during the Cooling Process

Figures5 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.Dislocation density distributions at room temperature (1, 3, and 5).
(a)
(b)
(c)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

Figures7 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
(10)σvon=(32SijSij)1/2.Residual stress distributions at room temperature (1, 3, and 5).
(a)
(b)
(c)We calculated the average residual stresses based on the plastic and elastic body models, as shown in Figure7. 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.Figure9 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.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)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.
(a)
(b)


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



---

*Source: 706923-2013-06-20.xml*