Computer modeling of crystal growth of silicon for solar cells

Abstract

A computer simulator with a global model of heat transfer during crystal growth of Si for solar cells is developed. The convective, conductive, and radiative heat transfers in the furnace are solved together in a coupled manner using the finite volume method. A three-dimensional (3D) global heat transfer model with 3D features is especially made suitable for any crystal growth, while the requirement for computer resources is kept permissible for engineering applications. A structured/unstructured combined mesh scheme is proposed to improve the efficiency and accuracy of the simulation. A dynamic model for the melt-crystal (mc) interface is developed to predict the phase interface behavior in a crystal growth process. Dynamic models for impurities and precipitates are also incorporated into the simulator.

Applications of the computer simulator to Czochralski (CZ) growth processes and directional solidification processes of Si crystals for solar cells are introduced. Some typical results, including the turbulent melt flow in a large-scale crucible of a CZ-Si process, the dynamic behaviors of the mc interface, and the transport and distributions of impurities and precipitates, such as oxygen, carbon, and SiC particles, are presented and discussed. The findings show the importance of computer modeling as an effective tool in the analysis and improvement of crystal growth processes and furnace designs for solar Si material.

Keywords

computer modeling / silicon / crystal growth / solar cells

Introduction

For several decades, the Czochralski (CZ) growth method has been the most widely used technique for growing high-quality silicon bulk crystals. Directional solidification (DS) is extensively employed in producing multi-crystalline silicon for solar cells because of its mass productivity and cost-effectiveness. Crystalline Si produced using these methods is the most commonly used material for solar cells [1]. The requirement for large and high quality wafers for the photovoltaic industry has made the design of the crystal growth furnace and the growth process a very challenging task. The thermal field, impurity distribution, melt-crystal (mc) interface, and melt flow motion have significant effects on the micro-defect formation and crystal quality in a crystal growth process. Thus, elucidating and controlling heat and mass transfer, impurity distribution, and the mc interface shape is important. With the development of modern computers and computation technology, numerical simulation has become an effective and essential tool for the design and optimization of a crystal growth process.

In the past two decades, extensive work has been devoted to the numerical analyses of bulk crystal growth processes using various models [2-11]. Most of these models can be generally divided into local models [2-5] and 2D global models [6-9]. However, the melt flow in a crucible and the thermal field within the growth furnace are usually three-dimensional (3D) under the influence of any asymmetric external fields or when using a square crucible in the DS. Therefore, a 3D global analysis of a growth system is necessary to gain a better understanding of growth phenomena and insight into the physics of growth processes. To this end, we develop a numerical simulator implemented with a 3D global heat transfer model, thereby making 3D global modeling feasible with moderate requirements for computer resources and computation time [10,11]. The developed simulator has been used successfully in quantitatively predicting the 3D features of a crystal growth. It has also been applied to a variety of crystal growth processes, such as solution growth for GaN [12,13], the PVT method for SiC [14,15], CZ growth, and DS processes for Si crystals [10,11,16-22]. To improve the efficiency and accuracy of the global simulation, we also develop a structured/unstructured combined mesh scheme. The large eddy simulation (LES) method is implemented in the simulator to simulate turbulence and the unstable melt flow in a large crucible. In the present study, the applications of the simulator to CZ growth and DS processes of Si crystals for solar cells are introduced.

Computation technologies of the numerical simulator for global simulation

A structured/unstructured combined multi-block mesh scheme for space discretization

Mesh quality is very important in numerical simulations. The use of a structured mesh enables the easy implementation of many algorithms of high computational efficiency. Meanwhile, the unstructured mesh is favorable for irregular geometries. In the current work, a structured/unstructured combined mesh scheme is developed to improve the efficiency and accuracy of global simulation. Figure 1 shows the configuration and computational mesh of a DS furnace. The argon flow field is discretized with an unstructured mesh. The other domains in the furnace are discretized with structure grids, as shown in Figs. 1(b) and (c). In the global simulation, all these block domains are solved in a coupled manner.

A 2D/3D combined multi-block mesh scheme for 3D global simulation

The typical configuration of a CZ furnace is illustrated in Fig. 2(a). All of the constituents of the furnaces are subdivided into a set of block regions. When 3D global modeling is to be conducted, a mixed 2D/3D finite volume scheme is used [10]. Following this scheme, the components in the core region of the furnaces (the 3D domain) are discretized in a three-dimensional manner, whereas the domains far from the core region (the 2D domain) are discretized in a two-dimensional manner. In Fig. 2, for example, the 3D domain includes the crystal, melt, crucible, and heater. The other regions in the furnace are included in the 2D domain. In the global simulation, the requirements for computation time and computer memory decrease considerably as all block domains are solved in a coupled manner.

LES for the melt flow with a large volume

The LES method is applied to simulate the unstable melt flow in a large crucible. The mass conservation, momentum, and energy equations for the melt flow in a crucible under the influence of a turbulence transport can be written as

(1)\nabla ·V=\mathrm{0},

(2)\rho V·\nabla V=-\nabla p+\nabla ·[{\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}(\nabla V+\nabla V^{\mathrm{T}})]-\rho g{\beta }_{\mathrm{T}}(T-T_{\mathrm{m}}),

(3)\rho c_{p}V·\nabla T=\nabla ·({\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}\nabla T),

where {\mu }_{\mathrm{e}\mathrm{f}\mathrm{f}}=\mu +{\mu }_{\mathrm{S}\mathrm{G}\mathrm{S}} and {\lambda }_{\mathrm{e}\mathrm{f}\mathrm{f}}=\lambda +{\mu }_{\mathrm{S}\mathrm{G}\mathrm{S}}/{\sigma }_{\mathrm{S}\mathrm{G}\mathrm{S}}, {\sigma }_{\mathrm{S}\mathrm{G}\mathrm{S}} are the turbulence Prandtle numbers, and {\mu }_{\mathrm{S}\mathrm{G}\mathrm{S}} is calculated using the dynamic Smagorinsky sub-grid scale model [23] as

(4){\mu }_{\mathrm{S}\mathrm{G}\mathrm{S}}=\mathrm{2}\rho C_{\mathrm{D}}{\overline{\mathrm{\Delta }}}^{\mathrm{2}}\left|\overline{S}\right|,

where \overline{\mathrm{\Delta }}={(\mathrm{V}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e})}^{\mathrm{1}/\mathrm{3}}, \left|S\right|=\sqrt{\mathrm{2}{\overline{S}}_{ij}{\overline{S}}_{ij}}, C_{\mathrm{D}}={\displaystyle \frac{L_{{\xi }_l{\xi }_m}{M}_{{\xi }_l{\xi }_m}}{M_{{\xi }_l{\xi }_m}{M}_{{\xi }_l{\xi }_m}}}, and L_{{\xi }_l{\xi }_m}=\tilde{\overline{u_{{\xi }_l}}}\tilde{\overline{u_{{\xi }_m}}}-\tilde{\overline{u_{{\xi }_l}}\overline{u_{{\xi }_m}}}, M_{{\xi }_l{\xi }_m}=\mathrm{2}{\overline{\mathrm{\Delta }}}_{\mathrm{C}}^{}\left[{\left({\widetilde{\overline{\mathrm{\Delta }}}}_{\mathrm{C}}/{\overline{\mathrm{\Delta }}}_{\mathrm{C}}\right)}^{\mathrm{2}}\left|{\tilde{\overline{S}}}_{\mathrm{C}}\right|{\tilde{\overline{S}}}_{{\xi }_l{\xi }_m}-\widetilde{\left|{\overline{S}}_{{\xi }_l{\xi }_m}\right|{\overline{S}}_{{\xi }_l{\xi }_m}}\right], and {\overline{S}}_{{\xi }_l{\xi }_m}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_{{\xi }_l}}{\partial {\xi }_m}}+{\displaystyle \frac{\partial u_{{\xi }_m}}{\partial {\xi }_l}}\right).

The global solution is obtained using an iterative procedure consisting of a set of local iterations for all block regions, with the calculation for radiative heat transfer in the furnace and a global conjugated iteration among them. The mc interface is calculated using a dynamic interface tracking method. Details on the models, including treatments of boundary conditions, can be found in Refs. [10,19].

Computer modeling of CZ-Si crystal growth processes

Dynamic behavior of the thermal field and mc interface of an industrial CZ-Si growth process

Simulations are carried out for an industrial CZ-Si crystal growth process. The diameter of the grown crystal is 200 mm. Figures 3(a)–(c) demonstrates the thermal fields in the core of the furnace at different growth stages. The dynamic behavior of the mc interface shape is presented in Fig. 3(d), which shows that the mc interface changes shape during the crystal growth process.

3D features of TMCZ-Si growth

Some results of the melt convection, thermal field, and mc interface shape are shown in Fig. 4 for a TMCZ configuration with a rotating crucible and crystal [16]. The transverse magnetic field is 0.1 T oriented in the x-direction. The rotation rates of the crystal and crucible are -30 and 5 r/min, respectively. Figure 4(a) shows the melt convection, thermal field, and mc interface profiles in symmetric planes x=0 (perpendicular to the magnetic field) and y=0 (parallel to the magnetic field). A strong vortex occupies plane x=0, whereas in plane y=0, the melt flows down to the crucible bottom with a complex pattern right under the mc interface. Accordingly, the temperature distributions in planes x=0 and y=0, which are dominated by the melt flow pattern, are distinctly different. However, even though the flow and thermal field of the melt are obviously three-dimensional, the mc interface profiles in the two orthogonal planes are almost symmetric. Figure 4(b) shows a 3D view of the mc interface shape and the temperature distribution on the melt top surface as well as on the mc interface.

Instability of the melt flow simulated using the LES method

Melt flow instability in the growth process for large diameter Si crystals is usually intense. To analyze the melt flow instability in such a crystal growth, global simulations are performed for an industrial CZ-Si crystal growth process, in which the diameter of the grown crystal is as large as 300 mm. The unstable melt flow is simulated using the LES method. Figure 5 displays the fluctuation fields of temperature and velocity. The figure shows that the instabilities of the melt flow under the mc interface and the melt flow at the melt surface are intensive.

Oxygen transport in the melt of a Si crystal growth with an EMCZ-TMF configuration

Some results of temperature and oxygen distributions in the molten silicon of an EMCZ-TMF growth configuration are presented in Fig. 6 [18]. The crystal and the crucible are rotated at -30 and 5 r/min, respectively. The transverse magnetic field is homogeneous with an intensity of 0.1 T in the x-direction. One electrode of the electric circuit is attached to the seed crystal while the other is attached to the melt surface. An electric current of 0.2 A is imposed to flow through the crystal and melt. Figure 6(a) shows the 3D temperature distribution in the melt. The 3D structure of the temperature distribution in the melt changes tremendously on account of the applied electric current. Figure 6(b) indicates the 3D distribution of oxygen concentration in the melt. Thin diffusion layers with large oxygen concentration gradients are observed close to the crucible walls and the melt free surface. The oxygen concentration, however, is fairly homogeneous inside the melt. The predicted values of oxygen concentration in the molten silicon are in good agreement with the experimental measurements of Hirata and Hoshikawa [24], demonstrating the accuracy of the simulation.

Computer modeling of DS processes for multi-crystalline Si ingot

3D features of the thermal field and solidification front surface in a DS-Si process

The configuration and dimensions of a small DS furnace for producing multi-crystalline silicon used for solar cells is shown in Fig. 1. The shape of the crucible is square. The argon flow rate is set to zero and the pressure is set to 0.8 atm. Some numerical results are presented in Figs. 7 and 8. Figure 7 shows the thermal fields in the furnace when half of the melt has been solidified. Figure 7(a) shows the global distribution of temperature in the furnace, the Si melt flow in the crucible, and the flow field of argon in the chamber. Figure 7(b) shows the 3D features of temperature distributions in the domains of the melt, crystalline ingot, crucibles, and pedestal. The mc interface corresponds to the isothermal surface of 1685 K. Figure 8 shows the mc interface shape, temperature distributions in the melt-crystal domain, and melt convective flow field. The interface is slightly concave to the crystalline ingot. The temperature gradient at the interface is approximately 8 K/cm in the solid and 1.3 K/cm in the melt. The melt flows toward the central axis at the melt surface with a maximal magnitude of an order of approximately 2 mm/s.

Impurity distribution and particle precipitation in a DS-Si process

Figure 9 shows some simulation results of impurity distribution in a grown Si ingot obtained from a DS process. Figure 9(a) displays the distribution of iron concentration in a solidified silicon ingot [20]. Areas with high iron concentration form at the top because of the segregation phenomenon, by which iron is segregated from the melt to the solid. Areas with high concentrations of iron also form close to the crucible walls. Such areas are formed by iron diffusion from the crucible walls. This observation is based on the small activation energy of iron diffusion in the solid phase of silicon.

A precipitation model of impurities is implemented in the simulator. Figure 9(b) shows the calculated distributions of substitutional carbon and generated SiC particles in a cross-plane of the solidified ingot in a slow cooling process [21]. Both the substitutional carbon and SiC particles are clustered in a narrow periphery-top region of the ingot.

Summary

A numerical simulator was developed with some effective computation technologies. Applications of this numerical simulator to CZ-Si crystal growth and DS processes for Si ingots for solar cells were introduced. Some typical results were presented, showing the importance and effectiveness of computer modeling and numerical simulation in analyzing and improving Si crystal growth processes.

References

[1] Li J F, Wang S C, Zhang M J, Ma L J. China Solar PV Report. Beijing: China Environmental Science Press, 2007 (in Chinese)

[2] Derby J J, Brown R A. Thermal-capillary analysis of Czochralski and liquid encapsulated Czochralski crystal growth I. Simulation. Journal of Crystal Growth, 1986, 74(3): 605–624

[3] Kobayashi N. Computational simulation of the melt flow during Czochralski growth. Journal of Crystal Growth, 1978, 43(3): 357–363

[4] Kakimoto K, Liu L J. Numerical study of the effects of cusp-shaped magnetic fields and thermal conductivity on the melt-crystal interface in CZ crystal growth. Crystal Research and Technology, 2003, 38(7, 8): 716–725

[5] Krauze A, Muiznieks A, Muhlbauer A, Wetzel T. Ammon W. Numerical 3D modelling of turbulent melt flow in large CZ system with horizontal DC magnetic field I. Flow structure analysis. Journal of Crystal Growth, 2004, 262(1-4): 157–167

[6] Atherton L J, Derby J J, Brown R A. Radiative heat exchange in Czochralski crystal growth. Journal of Crystal Growth, 1987, 84(1): 57–78

[7] Dupret F, Nicodeme P, Ryckmans Y, Wouters P, Crochet M J. Global modelling of heat transfer in crystal growth furnaces. International Journal of Heat and Mass Transfer, 1990, 33(9): 1849–1871

[8] Li M W, Li Y R, Imaishi N, Tsukada T. Global simulation of a silicon Czochralski furnace. Journal of Crystal Growth, 2002, 234(1): 32–46

[9] Kalaev V V, Evstratov I Yu. Makarov Yu N. Gas flow effect on global heat transport and melt convection in Czochralski silicon growth. Journal of Crystal Growth, 2003, 249(1, 2): 87–99

[10] Liu L J, Kakimoto K. Partly three-dimensional global modeling of a silicon Czochralski furnace I. Principles, formulation and implementation of the model. International Journal of Heat and Mass Transfer, 2005, 48(21, 22): 4481–4491

[11] Liu L J, Kakimoto K. Partly three-dimensional global modeling of a silicon Czochralski furnace II. Model application: Analysis of a silicon Czochralski furnace in a transverse magnetic field. International Journal of Heat and Mass Transfer, 2005, 48(21, 22): 4492–4497

[12] Kakimoto K, Liu L, Miyazawa H, Nakano S, Kashiwagi D, Chen X J, Kangawa Y. Numerical investigation of crystal growth process of bulk Si and nitrides-a review. Crystal Research and Technology, 2007, 42(12): 1185–1189

[13] Kashiwagi D, Gejo R, Kangawa Y, Liu L J, Kawamura F, Mori Y, Sasaki T, Kakimoto K. Global analysis of GaN growth using a solution technique. Journal of Crystal Growth, 2008, 310(7-9): 1790–1793

[14] Chen X J, Liu L J, Tezuka H, Usuki Y, Kakimoto K. Numerical investigation of induction heating and heat transfer in a SiC growth system. Crystal Research and Technology, 2007, 42(10): 971–975

[15] Chen X J, Liu L J, Tezuka H, Usuki Y, Kakimoto K. Optimization of the design of a crucible for a SiC sublimation growth system using a global model. Journal of Crystal Growth, 2008, 310(7-9): 1810–1814

[16] Liu L J, Kakimoto K. 3D global analysis of CZ-Si growth in a transverse magnetic field with rotating crucible and crystal. Crystal Research and Technology, 2005, 40(4, 5): 347–351

[17] Liu L J, Nakano S, Kakimoto K. An analysis of temperature distribution near the melt-crystal interface in silicon Czochralski growth with a transverse magnetic field. Journal of Crystal Growth, 2005, 282(1, 2): 49–59

[18] Liu L J, Nakano S, Kakimoto K. Investigation of oxygen distribution in electromagnetic CZ-Si melts with a transverse magnetic field using 3D global modeling. Journal of Crystal Growth, 2007, 299(1): 48–58

[19] Liu L J, Nakano S, Kakimoto K. Three-dimensional global modeling of a unidirectional solidification furnace with square crucibles. Journal of Crystal Growth, 2007, 303(1): 165–169

[20] Liu L J, Nakano S, Kakimoto K. Dynamic simulation of temperature and iron distributions in a casting process for crystalline silicon solar cells with a global model. Journal of Crystal Growth, 2006, 292(2): 515–518

[21] Liu L J, Nakano S, Kakimoto K. Carbon concentration and particle precipitation during directional solidification of multi-crystalline silicon for solar cells. Journal of Crystal Growth, 2008, 310(7-9): 2192–2197

[22] Miyazawa H, Liu L J, Hisamatsu S, Kakimoto K. Numerical analysis of influence of tilt of crucibles on interface shape and fields of temperature and velocity in a unidirectional solidification process. Journal of Crystal Growth, 2008, 310(6): 1034–1039

[23] Evstratov I Yu, Kalaev V V, ZhmakinA I,MakarovYu N,AbramovA G,IvanovN G,SmirnovE M,DornbergerE,VirbulisJ,TomzigE,von AmmonW. Modeling analysis of unsteady three-dimensional turbulent melt flow during Czochralski growth of Si crystals. Journal of Crystal Growth, 2001, 230(1, 2): 22–29

[24] Hirata H, Hoshikawa K. Oxygen incorporation and melt convection in CZ silicon crystal growth. Journal of the Japanese Association of Crystal Growth, 1988, 15(2): 207–216

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg