Nowadays there is an ever growing demand for a continuous development of the industrial silicon single crystal growth techniques. Complexity of the processes and high experimental costs is the motivating factor for advancing the mathematical modelling. We present our general approach for the modelling of industrial floating zone (FZ) process considering 3D features of molten zone shape and heat exchange there.

Considering the partly axisymmetric characteristics of the FZ system, axisymmetric phase boundary shapes, temperature field and radiation heat exchange (view factor method) in a quasi-stationary state are obtained by iteratively coupled calculations with program FZone [1]. The high frequency (≈3MHz) EM field is essentially asymmetric due to the asymmetry of the single turn inductor and therefore is modelled using 3D boundary element method [2]. Since the crystals are rotated, a reasonable assumption is to azimuthally average the obtained current density distributions for boundary conditions of the axisymmetric calculations. Illustrative calculation results as temperature field and shapes of phase boundaries for a typical 4” FZ crystal (process data from IKZ, Berlin) are depicted in Fig. 1.

Figure 1. Calculation results obtained with program FZone. Axisymmetric temperature field in feed rod, melt and crystal (left), isotherms in melt (right).

Nevertheless due to a pronounced asymmetry of the EM field for a more detailed considering of the molten zone, the influence of induced power density and EM pressure on silicon free melt surface cannot be regarded as axisymmetric. E.g., there are distinct maxima of induced power density on the free melt surface below the slits of inductor (Fig. 2) and therefore temperature field in the molten zone will be asymmetric. The maxima of electromagnetic pressure correspond to the maxima of the induced power density and will lead to asymmetry of the shape of the free melt surface.

Figure 2. Calculation results with program FZone. Induced current lines on silicon surfaces (left) and induced power density on the free melt surface (right).

We have developed a calculation algorithm for solving the equilibrium shape of free melt surface. The approach of our model is similar to that of the program Surface Evolver with energy functional minimization, i.e., finding equilibrium of gravity, surface tension, hydrostatic and EM forces which all are acting on the free surface. Triple point line (TPL) positions are fixed. Grid vertices are shifted in normal direction proportionally to the total force exerted on the vertex. Hydrostatic pressure p0 is chosen such that the average meniscus angle of the free surface at the crystal TPL is equal to the silicon growth angle. Since free melt surface shape affects the EM pressure distribution, solving both the EM and the free surface shape equations is iteratively coupled. Obtained calculation results for the shape of free melt surface, including comparison of azimuthal distributions of meniscus angle for 4” and 6” processes, is shown in Fig. 3.

Figure 3. EM pressure distribution (top left), isolines of vertical coordinate (top right) on the free surface of 4” crystal and azimuthal distribution of meniscus angle for 4” and 6” processes (bottom).

The asymmetry of induced power density creates asymmetric temperature field in the melt, which might affect the crystallization interface. The model for heat transfer in the molten zone is modelled as steady state heat diffusion, where temperature distribution is determined by Laplace equation with constant temperature boundary conditions on the crystallization and melting interfaces and given heat flux density on free melt surface. Direct boundary element method is used to solve the Laplace equation. Temperature values on the free surface and temperature normal derivatives on the crystallization interfaces are used as unknowns. The heat flux density on the free melt surface, which is used to formulate the boundary conditions, depends on the distributions of the EM power density, absorbed radiation heat flux density, and surface temperature. Due to the temperature dependence of the boundary conditions, the temperature equations are solved repeatedly, each time updating the boundary conditions until a steady solution is obtained. Results of the 3D heat exchange program are shown in Fig. 4.

Figure 4. View of the molten zone with boundary element mesh (top and centre), temperature distribution on free melt surface (bottom left) and temperature normal derivative on crystallization interface (bottom right).

The described mathematical models and respective programs have been incorporated into a calculation cycle with program FZone. The converged quasi-stationary axisymmetric solution obtained by FZone is used to generate a 3D mesh for iterative 3D EM field and free surface calculations. The calculated 3D free surface shape and EM power distribution on it are used to run 3D temperature calculations in the molten zone. The EM power distribution and absorbed radiated heat is used to formulate the boundary conditions for the free melt surface. The result of the 3D BEM heat transfer calculations is heat flux density distributions on the crystallization and melting interfaces. These distributions are azimuthally averaged and supplied to the FZone program for the use in 2D axisymmetric phase boundary calculations. The cycle is repeated until the necessary precision is reached.

Acknowledgements

The present work is carried out at the University of Latvia and has been supported by the European Regional Development Fund, project contract No. 2011/0002/2DP/2.1.1.1.0/10/APIA/VIAA/085.

References

[1] G. Ratnieks et al., Modelling of phase boundaries for large industrial FZ silicon crystal growth with the needle-eye technique, JCG, 255, 227-240 (2003)

[2] A. Muižnieks et al., Development of numerical calculation of electromagnetic fields in FZ crystal growth process, MHD, 46, 475-486 (2010)

• Legal notice:
  

  Copyright (c) Pielaszek Research, all rights reserved.
  The above materials, including auxiliary resources, are subject to Publisher's copyright and the Author(s) intellectual rights. Without limiting Author(s) rights under respective Copyright Transfer Agreement, no part of the above documents may be reproduced without the express written permission of Pielaszek Research, the Publisher. Express permission from the Author(s) is required to use the above materials for academic purposes, such as lectures or scientific presentations.
  In every case, proper references including Author(s) name(s) and URL of this webpage: https://science24.com/paper/28999 must be provided.

1. Hydrodynamical aspects of floating zone silicon crystal growth process
2. Numerical study of silicon crystal ridge growth