Article pubs.acs.org/IECR
Silicon Wafers for Solar Cells by Horizontal Ribbon Growth German A. Oliveros,† Ruochen Liu,† Seetharaman Sridhar,‡ and B. Erik Ydstie*,† †
Department of Chemical Engineering and ‡Department of Materials Science and Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, United States ABSTRACT: In this paper, we present a set of models that describe the principal components of the Horizontal Ribbon Growth processmainly, the interaction between fluid flow and heat transfer, the crystallization dynamics, and the effect of impurities on the morphology of the interface. Fluid-flow and heat-transfer models show the relationship between the pulling rate and the thickness of the silicon film. A crystallization model is developed to find the concentration distribution of impuritiesaluminum in this casein the melt and in the ribbon. We find that, because of low growth velocities, there is no formation of a soluteenriched boundary layer and that a 50-fold reduction of aluminum impurities can be expected. Finally, we use the Mullins− Sekerka stability theory to show that aluminum impurities at the proposed levels do not destabilize the interface upon applied perturbations. rate of crystallization.16 RGS decouples the direction of crystallization from the production rate and high pulling speeds can be obtained. However, the quality of the product is low, since the wafer tends to have very small crystalline domains. This is explained by the fact that the solid substrate upon which the melt is resting offers plenty of nucleation sites. Efficiencies greater than ∼10% have not been achieved on an industrial scale. The HRG process withdraws the silicon sheet horizontally, as in RGS, but the crystal grows vertically, as in the CZ process. In addition, the HRG process has the potential to produce monocrystalline crystals, since the ribbon rests on top of its melt; therefore, the complications arising in trying to obtain a monocrystalline crystal on a foreign substate are eliminated. It can therefore be expected that a high-quality wafer can be produced at high speed. Bleil13 verified the process concept by producing ice and germanium ribbons about a decade after Shockley15 had published the initial ideas. The paper by Bleil shows that the crystal thickness can be controlled by the temperature at the edge of the crucible and by manipulating the production rate. Inspired by Bleil’s experiments, Kudo14 proposed a similar design for HRG, using silicon as a substrate. A gas cooling system was introduced to convect heat away from the melt/crystal interface. Partial success was demonstrated in the investigation, and monocrystalline silicon sheets were produced at the rate of ∼20 cm/min for a limited period of time. However, the system was poorly optimized, difficult to control, and sensitive to disturbances and small changes in operating conditions. Different types of process instabilities were evident, such as uncontrolled dendritic growth and thickness variability. The silicon sheet produced by Kudo was uneven and unsuitable for solar applications. Nevertheless, the study was important, in that it was demonstrated, for the first time, that a silicon sheet could be produced by HRG. In a parallel study at the Jet
1. INTRODUCTION AND REVIEW OF LITERATURE Figure 1 shows that producing the silicon wafer accounts for ∼65% of the cost of a solar cell module.1,2 It is clear then that, in order to significantly reduce the cost of solar electricity, it is necessary to reduce the cost of wafer manufacture and silicon feedstock. In a previous paper, we studied how the fluidized-bed pyrolysis of silane can make less-expensive polysilicon feedstock than current batch technology.3 Other challenges and opportunities are reviewed by Gratzel,4 and a review of the supply chain from a process systems perspective is provided by Ranjan et al.2 The two-step Czochralski (CZ) process and the Bridgman (B) process contribute to more than 80% of the silicon substrate used in the solar cell industry.5 The CZ process produces a single crystalline cylindrical boule in Step 1. The Bridgman process produces a rectangular block of silicon. In Step 2, the silicon boule/block is sliced using wire saws to produce wafers that are ∼200 μm thick. The main advantage of the CZ and B processes is that they produce high-quality wafers on an industrial scale. The CZ process produces wafers that give a conversion efficiency of 18% or higher. Bridgman wafers typically have a conversion efficiency of ∼16%. The disadvantages of the CZ and B processes are many. Approximately 50% of the silicon raw material is lost during the sawing, doubling the cost, relative to direct wafering. Slow cooling to achieve single crystals (via the CZ process) or large crystalline domains (via the B process) give large energy costs. Finally, it is difficult to realize the economy of scale observed with batchwise operation. Scaleup of the CZ and B processes implies duplicating equipment, resulting in a nearly linear relationship between capital and operating expenses, and plant size.1 Many processes make silicon wafers continuously. Examples include Edge-Defined Film-Fed Crystal Growth (EFG),6−8 Dendritic Web (WEB),9,10 Ribbon Growth on Substrate (RGS),1,11,12 and Horizontal Ribbon Growth (HRG).13−15 Table 1 provides a comparison of the continuous methods, together with estimated production rates. The EFG and the RGS processes have been commercialized, but, so far, they have failed to produce high-quality wafers at low cost. The production rate of EFG is slow, since it is limited by the © 2013 American Chemical Society
Special Issue: Process Engineering of Energy Systems Received: Revised: Accepted: Published: 3239
July 14, 2012 January 14, 2013 January 15, 2013 January 15, 2013 dx.doi.org/10.1021/ie301857p | Ind. Eng. Chem. Res. 2013, 52, 3239−3246
Industrial & Engineering Chemistry Research
Article
Figure 1. Cost distribution for Czochralski-based silicon solar modules.
and surface-tension-driven (Marangoni) flow, in the presence of moving boundaries. Fairly complete process models have also been developed that include radiation cooling, heat transfer to the wall, and the transient effects.20,21 These studies form the foundation for the fluid flow models of the HRG process used in our study. We model the dynamics of the solidification and purification, with respect the aluminum impurities, by using the Stefan model in a moving (Lagrangean) frame of reference. The stability of the interface is investigated using Mullins−Sekerka theory.27
Table 1. Current Continuous Processes To Manufacture Silicon Ribbonsa
a
method
pull speed [cm/min]
width [cm]
throughput [cm2/min]
EFG WEB RGS vHGR
1.65 1−2 600−1000 10−20
8 × 12.5 5−8 12.5 5−8
5−16 165 7500−12500 160
Adapted from ref 1.
Propulsion Laboratory, Zoutendyk17 studied the influence of the pull velocity and the thermal boundary layer below the melt/ crystal interface. HRG was apparently abandoned as a way to produce solar cells once the oil crisis of the 1970s subsided and the interest in low-cost, ground-based photovoltaics dissipated. The experimental studies by Bleil and Kudo show that it is necessary to develop a better understanding for the static and dynamic characteristics of silicon crystallization in HRG. Recent studies2,18 show that is necessary to model the interconnections among solid and fluid flow, heat transfer, and crystal growth since instabilities originate at the solidification interface, as shown in Kudo’s experimental program.14 It is also necessary to study how impurities are removed and how they affect process stability to demonstrate the purification aspect of the process claimed by Shockley. Silicon is one of the most-studied elements, because of its importance as a semiconductor material. It is not surprising then that silicon crystal growth has been investigated quite extensively.19 Very comprehensive model systems have been developed for the CZ process, using discretization in finite elements. These models include buoyancy-driven (Boussinesq)
2. FLUID FLOW AND HEAT TRANSFER In this section, we develop a two-dimensional (2D) steady-state fluid flow and heat-transfer model of HRG to estimate the heat flow from the melt to the cooling medium through the solid silicon film. Figure 2 shows a version of Shockley’s15 conceptual design of a silicon wafer process. Small pellets of silicon melt in Zone I. The molten silicon then flows to the stabilization zone (Zone II), where steady flow and temperature are assumed. In Zone III, the silicon cools and a sheet of silicon forms on top of the melt. The cooling and production rates control the sheet thickness as it continuously exits the process on the right. The sheet is then cut and processed further to form solar cells. At steady state, crystallization mostly occurs in the vertical direction, as illustrated in models and experiments by Bleil and Kudo. Since we are mainly interested in the relationship between pulling rate and thickness of the ribbon at the steady state, the effects of thermal undercooling and interface growth and stability will be studied in detail in a separate model. This is done to reduce the computational complexity of both models and, at the same time, be able to gather as much information from the system as
Figure 2. Conceptual design of the silicon wafering process. 3240
dx.doi.org/10.1021/ie301857p | Ind. Eng. Chem. Res. 2013, 52, 3239−3246
Industrial & Engineering Chemistry Research
Article
possible. A comprehensive model of the HRG is found in the work by Daggolu, Derby, and co-workers.18 The HRG process is somewhat similar to the Pilkington float process used to form flat glass.22 The main difference is that a float substrate, such as tin, is not needed, since solid silicon floats on its liquid, just like ice floats on water. One major practical challenge with the process is that it is very difficult to perform experiments with molten silicon. The operating temperature is high and silicon is very reactive, so it is difficult to design and maintain equipment to control the process and ensure that silicon retains the purity specifications needed to make high-efficiency solar cells. Some other challenges are outlined by Kudo14 and Zoutendyk.17 Figure 3 shows the geometry of the system that we have modeled using the COMSOL Multiphysics modeling environ-
∇·(k∇T (x , y)) = ρCpu(x , y) ·∇T (x , y)
Heat is released at the surface through the radiation boundary condition 4 k∇T (x , S) = εσ(Tamb − T (x , S)4 )
(5)
where the emissivity (ε) accounts for deviation from blackbody radiation and S represents the surface. The other boundary conditions include no slip at the walls and continuity of stress at the surface of the melt. The parameters used to model this system, and the rest of the systems in the paper, are listed in Table 2. Table 2. Material Properties and Parameters Used in All of the Mathematical Models
Figure 3. Simplified geometry of Zone II and Zone III of a silicon wafering system.
ment.23 It is assumed that silicon behaves as a Newtonian incompressible fluid.18,20,21 The continuity equation in two dimensions then reads ∇·u(x , y) = 0
(4)
(1)
where u is the velocity field. The momentum balance is then represented by the Navier−Stokes equations with a body force F to account for the phase change, so that we have24
parameter
symbol
value
melting temperature averaged liquid temperature gradient average crystal temperature gradient initial interface velocity conductivity of solid silicon conductiviy of liquid silicon density of solid silicon density of liquid silicon viscosity of silicon heat capacity of solid silicon thermal diffusivity of liquid silicon thermal diffusivity of solid silicon mass diffusivity Al−Si emissivity Stefan−Boltzmann constant latent heat of fusion segregation coefficient of Al−Si capillary constant slope of the liquidus line
Tmelt G G′ V ks kl ρs ρl μ Cp αl αs Dl ϵ σ ΔH kseg Γ m
1687 K 4.48 × 10−3 K m−1 4.49 × 10−3 K m−1 0.6 × 10−6 m s−1 18 W m−1 K−1 58 W m−1 K−1 2293 kg m−3 2570 kg m−3 0.58 × 10−3 Pa s−1 1040 J kg−1 K−1 2.33 × 10−5 m2 s−1 7.54 × 10−6 m2 s−1 7 × 10−8 m2 s−1 0.64 5.67 × 10−8 W m−2 K4 1.79 × 106 J kg−1 2.83 × 10−3 8.99 × 10−11 m −9.5 × 10−5 K ppm−1
ρ(u(x , y) ·∇u(x , y)) = ∇·[−p(x , y)I + η(∇u(x , y) + (∇u(x , y))T )] + F
The simulation shows that silicon solidifies on top of the molten substrate at steady state. Enough heat is released for the surface to reach the melting temperature (Figure 4) to form a sheet of ∼300 μm with a pulling rate of 16 cm/min. The sensitivity plot in Figure 5 shows that thickness decreases as the pulling rate is increased, as shown in the experiments by Bleil and Kudo.
(2)
The parameter η represents the viscosity of silicon and the body force F accounts for the “phase distinction”. It is defined so that25 ⎡ (1 − B)2 ⎤ ⎥Au F=⎢ 3 ⎣ B +q ⎦
(3)
In this expression, B represents the liquid fraction. The number q (∼0.001) prevents small divisors. The function F vanishes when B is equal to one (only fluid) and becomes large as B tends toward zero. Large F forces the velocity in the solid domain to be nearly uniform, because of the very high body force, and this mimics the behavior of a solid. This model does not provide a sharp interface between liquid and solid silicon. Instead, we get what is referred to as the “mushy zone”. By doing this, we avoid the issue of calculating the exact shape of the moving boundary, while still being able to caputure the desired operating parameters (thickness and pulling rate). In the next section, we develop a more-detailed model to locate the position of the freezing front, and we study the effect of impurities on the stability of the crystallization front. Molten silicon is fed at 1500 °C from the left boundary in Figure 3. The temperature distribution in the liquid is given by the energy transport equation,
Figure 4. Temperature profile of the process system and the velocity field of the molten silicon. 3241
dx.doi.org/10.1021/ie301857p | Ind. Eng. Chem. Res. 2013, 52, 3239−3246
Industrial & Engineering Chemistry Research
Article
The analysis below calculates the extent of purification achieved by HRG. The formation of the impurity layer at the liquid/solid interface can be detrimental for the process, since it can cause steep concentration gradients that induce interfacial instabilities (the growth of dendrites), as observed experimentally by Kudo.14 The model is developed from a Lagrangean perspective: as crystallization progresses, the growing interface “sweeps” the liquid domain and can be imagined as an observer sitting on top of the moving ribbon. The molten bath is assumed to be still, with respect to the film. We start the model development using silicon with aluminum with a bulk concentration of 50 ppm. This corresponds to the concentration that has been achieved in a metallurgical process for co-producing silicon and a silicon-rich aluminum alloy from kaolin. The model assumes that transport at the surface is dominated by diffusion. This assumption is justified since the fluid flow model developed above shows that only a thin layer of the melt is affected. In the liquid domain, the transient heat transfer and mass diffusion equations are given by the following partial differential equations (PDEs):
Figure 5. Sensitivity analysis of the computational fluid dynamics (CFD) model.
3. INTERFACE DYNAMICS One unique advantage of HRG is that it can be used as a purification processas any other crystallization process can for elements that have large difference in the solubility between the melt and the crystal, while, at the same time, avoiding contamination from carbon coming from a substrate or a die, such as in the RGS and EFG processes. The impurities diffuse preferentially to the liquid phase, leading to the formation of a solute-enriched boundary layer close to the crystallization front (see Figure 6). As solidification proceeds, solute redistributes
∂Tl ∂ 2T = αl 2l ∂t ∂y
(7)
∂Cl ∂ 2C = Dl 2 l ∂t ∂y
(8)
Impurity diffusion is negligible in the solid. Thus, the PDE for impurities in the solid is given by
∂Ts ∂ 2T = αs 2s ∂t ∂y Cs = ksegCl
(9) (10)
The energy balance dictates that the difference between the heat conducted from the liquid bath to the interface and the heat conducted from the film to the surface of the system balances the heat of crystallization released at the interface. Since the diffusion coefficient in the solid is much smaller than in the liquid (Ds < < Dl), it is assumed that solute diffusion is negligible and, therefore, the solute concentration in the solid is given by eq 13, where kseg is the equilibrium partition coefficient. An additional relationship between the thermal and concentration fields is required, since the melting point of the aluminum−silicon alloy is dependent on the amount of aluminum impurity. These three constraints are expressed as Energy Balance:
Figure 6. Solute is rejected at the interface due to the large difference in impurity solubility between the two phases.
among the two phases, depending on the magnitude of segregation coefficient (kseg) and the diffusivities of the impurities. The velocity of the crystallization front determines the extent of the diffusion boundary layer at the interface.26 Diffusion in the melt is uniform throughout the sample, and equilibrium-type behavior can be observed at low velocities. The phase diagram then determines the behavior at every point in the system. The solute builds up close to the interface at high velocities and a boundary layer develops. The boundary layer thickness is dependent on the interface velocity, and the liquid diffusivity is given by the equation D δ= l (6) V Large interface velocities “trap” the solute close to the interface, inducing depression of the melting point. (See Figure 6.) This undercooling provides the driving force for the development of morphological instabilities, as we explain in the next section.
ks
∂Ts ∂T ∂y − kl l − ρs ΔH =0 ∂y ∂y ∂t
(11)
Mass Balance: −Dl
∂Cl ∂y − (Cl − Csint) = 0 ∂y ∂t int
(12)
Melting Point for Binary Mixture: Tint = Tmelt + mCl
(13)
The term ∂y/∂t in the energy and mass conservation equations represents the velocity of the crystallization front; m defines the slope of the liquidus line in the aluminum−silicon phase diagram. Equation 13 provides an estimate for freezing-point depression, 3242
dx.doi.org/10.1021/ie301857p | Ind. Eng. Chem. Res. 2013, 52, 3239−3246
Industrial & Engineering Chemistry Research
Article
due to the accumulation of impurities close to the liquid/solid interface. In the following, we assume that the slope of the liquidus line is constant. The boundary condition at the surface of the silicon melt remains the same as in the fluid dynamics model. k
∂T = εσ(T 4 − Tamb 4) ∂y
(14)
Once the surface temperature reaches the solidification temperature at time tcrystallization (Figure 7), which is dependent
Figure 8. Evolution of the temperature profile in both phases close to the crystallization front.
Figure 10 shows the concentration profile of aluminum in molten silicon close to the crystallization front. The velocity of the crystallization front is shown in Figure 11. As the velocity of the crystallization front converges (dy/dt ≈ 0), the length of the solute-enriched boundary layer increases and the concentration profile “flattens”. For the proposed operating parameters, the crystallization rate varies between 0.9 and 0.1 μm per second, whereas the pull rate is 16 cm/min (2.6 × 103 μm/s). This shows that the rate of crystallization is much slower than the pull rate and provides an indication that a high-quality crystalline sheet can be produced. Figure 12 shows the solute distribution across the thin film. Liquid concentration from 50 ppm aluminum in the melt gives a 380-μm-thick film with an aluminum impurity of 0. In this expression, V is the velocity of the interface; G′ is the temperature gradient of the crystal at the interface, (∂Ts/∂y), multiplied by an average conductivity, k̅ = 0.5 × (ks + kl); likewise G is the averaged temperature gradient of the melt at the interface; Gc is the concentration gradient at the interface; and p = 1 − kseg. Finally, 3244
dx.doi.org/10.1021/ie301857p | Ind. Eng. Chem. Res. 2013, 52, 3239−3246
Industrial & Engineering Chemistry Research
Article
in the steep decrease of the perturbation function for values of ω less than unity. As the frequency increases, dominance of thermal gradients over surface energy increases. From the analysis above, we expect that the HRG process produces a flat front and a homogeneous crystal with an aluminum impurity of
