- Email: [email protected]
Numerical studies of flow patterns during Czochralski growth of square-shaped Si crystals
Journal of Crystal Growth 318 (2011) 244–248
Contents lists available at ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro
Numerical studies of flow patterns during Czochralski growth of square-shaped Si crystals W. Miller , Ch. Frank-Rotsch, P. Rudolph Leibniz Institute for Crystal Growth (IKZ), Max-Born-Str. 2, 12489 Berlin, Germany
a r t i c l e i n f o
a b s t r a c t
Available online 26 October 2010
We have performed 3D calculations of the melt flow during Czochralski growth of nearly square-shaped silicon crystals by means of Ansys-cfx. Beforehand, the temperature field in the furnace has been computed by axisymmetric computations using CrysMAS for providing thermal boundary conditions for the local simulation of the melt flow. CrysMAS has been also used for the computations of Lorentz force
Keywords: A1. Computer simulation A1. Fluid flows A2. Magnetic field assisted Czochralski method
s
g densities, produced by a travelling magnetic field using the KRISTMAG approach. & 2010 Elsevier B.V. All rights reserved.
1. Introduction In order to reduce material loss during cutting Si single crystals for photovoltaic wafers a quadratic shape would be favourable. It has been shown that the growth of Si single crystals with nearly quadratic cross section is possible by the float zone using a special HF inductor [1]. Recently, a technique has been developed to grow Si single crystals with nearly quadratic cross section by the Czochralski method using travelling magnetic fields (TMFs). The key point is a low radial but stable temperature gradient being the precondition for the appearance of four large rectangularly arranged {1 1 0} facets parallel to the pulling axis [0 0 1]. Experimentally, this has been achieved by applying an appropriate TMF with the heater magnet g s project [2]. The module (HMM) as developed in the KRISTMAG details of the experiments are presented in [3]. In this paper we will analyse the results of 3D calculations of the melt flow using the software Ansys-cfx. The thermal boundary conditions were obtained from a global axisymmetric calculation with CrysMAS [4]. Due to the high Reynolds numbers for such Czochralski system the flow computation is still challenging (see e.g. [5]). We use the large eddy simulations (LES) as it has become common for computing the flow in Czochralski growth of silicon [6,7].
2. Numerical methods We use the commercial software package Ansys-cfx for the 3D calculation of the melt flow. Ansys-cfx is a finite volume code with the algebraic multigrid technique (www.cadfamily.com/down load/CAE/CFX/xthry.pdf). Corresponding author. Tel.: +49 30 63923074; fax: + 49 30 63923003.
E-mail address: [email protected] (W. Miller). 0022-0248/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2010.10.098
Because the crystal is the only 3D object we let the crystal at rest and use a rotational frame with the angular velocity oX ¼ 2pOX , where OX is the rotation rate of the crystal. The angular velocity enters the Corriolis and centrifugal force. The boundary condition for the velocity at the crucible is set as vbc ¼ 2prðOc þ OX Þ, where OX is the rotation rate of the crucible.The procedure has been verified by computations for cylindrical crystals, where we computed the flow field also in the usual way. In the cases considered the Reynolds numbers are large and thus very small flow structures can be expected. Therefore, we used the large eddy Smagorinsky model. Here, the viscosity in the Navier– Stokes equation is given by mNS ¼ m þ msgs , where m is the viscosity of the melt (see Table 1) and 2 qffiffiffiffiffiffiffiffiffiffiffiffi 2Sij Sij ð1Þ msgs ¼ min fm Cs V 1=3 ,lmax is the turbulent viscosity. Sij is the shear stress of the velocity field, V the volume of the cell and Cs the Smagorinsky constant, which set to Cs ¼0.1. The factor fm takes the van Driest damping at the walls into account: e=AÞ3 , fm ¼ 1exp½ðy
ð2Þ
e ¼ rðywall u e : local veloe Þ=m (u with the normalized wall distance y city) and the damping factor A ¼25. lmax presents a cutoff in the case of volumes with large aspect ratios near the wall, i.e. with ywall 5 V 1=3 : lmax ¼ kywall . k was chosen to be 0.4. For all 3D calculations we used a hexahedral mesh with 931,000 volumes. A finer mesh was used at the boundaries up to a resolution of 0.12 mm at the first volume. Time steps have been chosen such that the root-mean-square Courant numbers were smaller than 1.
3. Characteristic numbers for convection The strength of the buoyancy convection is described by the Grashof number Gr ¼ g ar2 m2 l30 DT, where we chose the crucible
W. Miller et al. / Journal of Crystal Growth 318 (2011) 244–248
radius Rc as the characteristic length l0 and the maximum temperature in the melt (39 K) as DT. The other parametersp are ffiffiffiffiffiffigiven in Table 1. For analogy reasons one can identify Re ¼ Gr, where Re ¼ l0 u0 rm1 is the Reynolds number with a characteristic velocity u0. Thus, we obtain Rebuoy ¼ 1:8 104 :
ð3Þ
The Reynolds numbers for crucible and crystal rotation are defined by Rec ¼ 2pOc R2c rm1 and Rec ¼ 2pOX R2c rm1 , respectively. The rotation rates were OX ¼ 8:5 rpm and Oc ¼ 5 rpm for the crystal and crucible, respectively. With these values we obtain ReX ¼ 1:7 104 ,
Rec ¼ 9:8 103 :
ð4Þ
In the same way as for the buoyancy forces we can define a Reynolds number for the Lorentz forces as ReL ¼
maxðFL ÞR3c
rm2
:
ð5Þ
The maximum Lorentz force density max (FL) is about 2000 N m 3 (see below) and thus we get ReL ¼ 6 104 ,
ð6Þ
Table 1 Physical parameter used for the computations. Melting point temperature Density Specific heat Thermal conductivity Dynamic viscosity Thermal expansion factor Gravitational acceleration
Tm
1685 K 2420 kg m 3 1000 J kg 1 K 1 60 W m 1 K 1 7 10 4 Pa s 1.41 10 4 K 1 9.81 m s 2
r cP
l
m a g
245
which is significantly larger than any of the other Reynolds numbers. From u0 ¼ Rem=ðRc rÞ we can estimate the velocity range in a vertical plane of the melt, yielding u0 ¼ 230 mm/s. Indeed, we observed maximum velocities in this range for the run with TMF.
4. Results 4.1. Temperature field The experimental set-up with a heater magnet module (HMM) around the crucible consisting of three side coils and two bottom spirals is described in detail in [3]. For the numerical simulation we consider a growth situation, where the height of the melt was 41 mm (see Fig. 2). First of all, we computed the thermal field in the entire furnace using the axisymmetric code CrysMAS [4]. Applying the heater power as in the experiment we observed a very small radial temperature gradient. We like to note that this calculation did not include melt convection. The low temperature difference between crystal and crucible of about 9 K and in the entire melt of about 39 K is a result of the high thermal isolation. For the axisymmetric calculation we chose the maximal diameter of the crystal (see Fig. 1), i.e. 94 mm. The temperature data of CrysMAS are used as a boundary condition for the free melt surface in the 3D calculation. Because the temperature is almost the melting point temperature around the crystal it is not a large error when using the data of the axisymmetric simulation in 3D though the crystal is now not axisymmetric any longer. The temperature profile along the free interface of the melt can be seen in detail in Fig. 3. The crystal/ melt interface was computed to be slightly convex near to a flat interface. In the 3D calculations we use a completely flat crystal/ melt interface and neglect the wetting of crystal and crucible for simplifying the geometry and the grid generation. The crystal shape is taken as measured after the experiment (see Fig. 1): 94 mm diagonal diameter, 58 mm facet length. 4.2. Velocity field without TMF
94 mm
38 mm
152 mm
20 mm
crystal melt free surface
Fig. 1. Geometry of crystal and top of crucible as used for the computations. The data correspond to those of the crystal grown.
We start with 1/10 of the real Reynolds number and then increase the Reynolds numbers twice until the real values. For the runs with the real Reynolds numbers we use a time step of 0.01 s for the computation with the full Reynolds numbers. We observe the typical convection roll going up at the crucible’s side walls. In Fig. 7 (right) the averaged velocity field in a vertical cut is shown. Such a velocity field causes a heat transport from the sides of the crucible towards the crystal. From Fig. 6 (middle) it can be clearly seen that the radial temperature gradient near the free surface is increased. Because we fixed the temperature at the free surface this effect cannot be seen there. Despite this restriction, it is clear that maintaining a stable low radial thermal gradient would be not possible by crucible rotation. The other effect of the heat transport by the convection roll is the low temperature at the bottom of the crucible causing the danger of freezing (see Fig. 6).
+
+
+
+
39 mm
47 mm
74 mm
z = 30 mm
20 mm
z = 41 mm
z = 0 mm Fig. 2. Side view of geometry with locations of monitor points. The monitor point 2 is at the same distance from the axis as the {1 1 0} facet of the crystal, and point 3 is at the same distance as the diagonal section of the crystal. The monitor points are in a plane normal to the {1 1 0} facets.
-190
Fr / Nm−3
-1930
-1550
z
-1160
W. Miller et al. / Journal of Crystal Growth 318 (2011) 244–248
-770
246
r Fig. 3. Temperatures at crystal/melt interface and free melt interface as computed by CrysMAS. In the vicinity of the crystal the temperature is almost the melting point temperature, i.e. the radial gradient is nearly zero.
Fig. 5. Vectors of Lorentz force distribution and radial component Fr of the Lorentz force density.
1724
1715
1705
1695
1685
T [K]
Fig. 6. Comparison of the time-averaged (7 s) temperature field in a vertical plane for heat conduction only (top), buoyancy flow and rotation (middle), and buoyancy flow, rotation and TMF (bottom). Vertical plane is the same as in Fig. 7. Fig. 4. Instantaneous velocities on the free melt surface for the case without TMF. Thick arrows indicate the direction of crucible and crystal rotation.
At the free melt surface we can observe a recirculation of the flow behind the rotating edge of the crystal (see a snapshot in Fig. 4). The radially inwards and outwards directed flow is changing with time whereas the recirculation remains. The eddy viscosity is with a maximum value of msgs ¼ 1 104 Pa s significantly smaller than those of the melt, stating that the effects of numerically unresolved turbulent mixing is small. 4.3. Velocity field with TMF The Lorentz forces have been computed with CrysMAS and transferred to Ansys-cfx. The Lorentz forces are quite large with a maximum of about 2000 N m 3. Most important for the success of growing a Si crystal with {1 1 0} facets are the large radial components Fr in the lower side part of the crucible and the vertical components
at the bottom Fz directing upwards, i.e. opposite to the gravity (see Fig. 5). Such a Lorentz force distribution causes a convection roll at the crucible moving downwards at the side walls and upwards in the region of the crystal. The inversion of the convection roll by a TMF was demonstrated by Klein et al. in numerical simulations [8]. Here, we use much higher Lorentz forces yielding a Reynolds number more than four times larger than that of the buoyancy force and than that of the crucible rotation (see Section 3.2). Therefore, it can be expected that the flow induced by TMF is dominant. Indeed, the convection roll is completely inverted (see left side of Fig. 7) with a maximum time-averaged velocity of 115 mm/s compared to 15 mm/s for the case without TMF. Due to the higher velocities we had to use smaller time step (0.005 s) than in the case without TMF. Because the azimuthal velocities are now small compared to the vertical and horizontal velocities we did not observe a specific flow pattern due to the nearly quadratic shape of the crystal as in the case without TMF. The eddy viscosity is now up to msgs ¼ 4 103 Pa s in
W. Miller et al. / Journal of Crystal Growth 318 (2011) 244–248
247
Fig. 7. Time-averaged velocities (7 s) in a vertical plane with TMF (left half) and without TMF (right half). The lengths of the arrows correspond to the velocity but the arrow length is three times larger in the right half than in the left one. Maximum time-averaged velocity with TMF was 115 mm/s and without TMF 15 mm/s.
TMF 1690
TMF
1690
1688 1688
T [K]
1686 0
1
2
3
4
5
6
0
7
1
2
3
4
2
3
4
5
6
7
1706 1692 TMF
1702 1698
1690
1694 TMF 1690
1688 0
1
2
3
4
5
6
7 t [s]
0
1
5
6
7
Fig. 8. Temperature oscillations at the monitor points 1, 2, 3 and 4 as defined in Fig. 2. Case without TMF is red and case with TMF is blue as marked. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
certain regions, which is about six times the viscosity of the melt. A detailed analysis of the oscillating flow field and the distribution of the eddy viscosity will be published elsewhere. 4.4. Comparison of case without and with TMF We will compare the results without and with TMF in more detail. Firstly, we look on the time-averaged velocities in Fig. 7. In the case with TMF (left) we have a large convection roll fitting into the melt region below the free surface with a maximum velocity (time-averaged) of 115 mm/s. The flow beneath crystal is in the same direction as caused by crystal rotation and no shear flow occurs as in the case without TMF (right). In total, there is only a weak convection below the crystal. Without TMF we have the typical convection roll going upwards at the rim of the melt plus a counter vortex at the top outer part. As a result of the completely different velocity fields we obtain different temperature fields (see Fig. 6). The thermal gradient radial direction is much higher without TMF than with TMF, whereas the vertical gradient in the vicinity of the crystal/melt interface is higher in the case with TMF. The latter will stabilize a
flat crystal/melt interface. There is also no danger of freezing at the bottom in the case with TMF. The numerical calculations explain the experimental experience that the growth of a facetted crystal with a flat interface is possible using such a TMF. The time dependence at the four monitor points defined in Fig. 2 is given in Fig. 8 for both the case without TMF and that with TMF (marked). Oscillations with frequencies higher than 10 Hz have been suppressed by time-averaging. The following can be noticed:
The temperatures in the case with TMF are nearly the same at all monitor points.
The frequency is higher in the case with TMF, most pronounced at monitor point 1.
The oscillation amplitude is smaller in the case with TMF. 5. Conclusions We have performed 3D simulations of the melt flow with and without TMF. The TMF as used in the experiments causes a strong vortex with a heat transport from the crucible’s bottom to the crystal. Such a heat transport stabilizes the low temperature gradient along
248
W. Miller et al. / Journal of Crystal Growth 318 (2011) 244–248
the free surface and the stable growth of the silicon crystal with {1 1 0} facets becomes possible. The TMF induced vortex is in opposite direction to the common vortex caused by buoyancy convection and enhanced by crucible rotation. The results of the numerical calculations show that this rotation-driven vortex leads to a significant increase of the radial gradient near the free surface, which would inhibit the formation of {1 1 0} facets because of the decreased undercooled regions along the {1 1 0} places (see [3]). References [1] A. Muizˇnieks, A. Rudevics, K. Lacis, H. Riemann, A. Ludge, F.W. Schulze, B. Nacke, Square-shaped silicon crystal rod growth by the fz method with special 3D shaped HF inductors, Magnetohydrodynamics 43 (2) (2007) 269–282. [2] P. Rudolph, Travelling magnetic fields applied to bulk crystal growth from the melt: the step from basic research to industrial scale, J. Crystal Growth 310 (2008) 1298–1306.
[3] P. Rudolph, M. Czupalla, B. Lux, F. Kirscht, C. Frank-Rotsch, W. Miller, M. Albrecht, The use of heater-magnet module for Czochralski growth of PV silicon crystals with quadratic cross section, J. Crystal Growth, this issue, doi:10.1016/j.jcrysgro.2010.10.070. [4] J. Fainberg, D. Vizman, J. Friedrich, G. Mueller, A new hybrid method for the global modeling of convection in CZ crystal growth configurations, J. Crystal Growth 303 (2007) 124–134. [5] W. Miller, Numerical simulations of bulk crystal growth on different scales: silicon and GeSi, Phys. Status Solidi B 257 (4) (2010) 855–869. [6] A. Muiznieks, A. Krauze, B. Nacke, Convective phenomena in large melts including magnetic fields, J. Crystal Growth 303 (2007) 211–220. [7] I.Y. Evstratov, V.V. Kalaev, A.I. Zhmakin, Y.N. Makarov, A.G. Abramov, N.G. Ivanov, A.B. Korsakov, E.M. Smirnov, E. Dornberger, J. Virbulis, E. Tomzig, W. von Ammon, Numerical study of 3D unsteady melt convection during industrial-scale CZ Si-crystal growth, J. Crystal Growth 237–239 (3) (2002) 1757–1761. [8] O. Klein, C. Lechner, P.-E. Druet, P. Philip, J. Sprekels, C. Frank-Rotsch, F.-M. Kießling, W. Miller, U. Rehse, P. Rudolph, Numerical simulation of Czochralski crystal growth under the influence of a traveling magnetic field generated by an internal heater-magnet module (HMM), J. Crystal Growth 310 (2008) 1523–1532.
Contents lists available at ScienceDirect
Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro
Numerical studies of flow patterns during Czochralski growth of square-shaped Si crystals W. Miller , Ch. Frank-Rotsch, P. Rudolph Leibniz Institute for Crystal Growth (IKZ), Max-Born-Str. 2, 12489 Berlin, Germany
a r t i c l e i n f o
a b s t r a c t
Available online 26 October 2010
We have performed 3D calculations of the melt flow during Czochralski growth of nearly square-shaped silicon crystals by means of Ansys-cfx. Beforehand, the temperature field in the furnace has been computed by axisymmetric computations using CrysMAS for providing thermal boundary conditions for the local simulation of the melt flow. CrysMAS has been also used for the computations of Lorentz force
Keywords: A1. Computer simulation A1. Fluid flows A2. Magnetic field assisted Czochralski method
s
g densities, produced by a travelling magnetic field using the KRISTMAG approach. & 2010 Elsevier B.V. All rights reserved.
1. Introduction In order to reduce material loss during cutting Si single crystals for photovoltaic wafers a quadratic shape would be favourable. It has been shown that the growth of Si single crystals with nearly quadratic cross section is possible by the float zone using a special HF inductor [1]. Recently, a technique has been developed to grow Si single crystals with nearly quadratic cross section by the Czochralski method using travelling magnetic fields (TMFs). The key point is a low radial but stable temperature gradient being the precondition for the appearance of four large rectangularly arranged {1 1 0} facets parallel to the pulling axis [0 0 1]. Experimentally, this has been achieved by applying an appropriate TMF with the heater magnet g s project [2]. The module (HMM) as developed in the KRISTMAG details of the experiments are presented in [3]. In this paper we will analyse the results of 3D calculations of the melt flow using the software Ansys-cfx. The thermal boundary conditions were obtained from a global axisymmetric calculation with CrysMAS [4]. Due to the high Reynolds numbers for such Czochralski system the flow computation is still challenging (see e.g. [5]). We use the large eddy simulations (LES) as it has become common for computing the flow in Czochralski growth of silicon [6,7].
2. Numerical methods We use the commercial software package Ansys-cfx for the 3D calculation of the melt flow. Ansys-cfx is a finite volume code with the algebraic multigrid technique (www.cadfamily.com/down load/CAE/CFX/xthry.pdf). Corresponding author. Tel.: +49 30 63923074; fax: + 49 30 63923003.
E-mail address: [email protected] (W. Miller). 0022-0248/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2010.10.098
Because the crystal is the only 3D object we let the crystal at rest and use a rotational frame with the angular velocity oX ¼ 2pOX , where OX is the rotation rate of the crystal. The angular velocity enters the Corriolis and centrifugal force. The boundary condition for the velocity at the crucible is set as vbc ¼ 2prðOc þ OX Þ, where OX is the rotation rate of the crucible.The procedure has been verified by computations for cylindrical crystals, where we computed the flow field also in the usual way. In the cases considered the Reynolds numbers are large and thus very small flow structures can be expected. Therefore, we used the large eddy Smagorinsky model. Here, the viscosity in the Navier– Stokes equation is given by mNS ¼ m þ msgs , where m is the viscosity of the melt (see Table 1) and 2 qffiffiffiffiffiffiffiffiffiffiffiffi 2Sij Sij ð1Þ msgs ¼ min fm Cs V 1=3 ,lmax is the turbulent viscosity. Sij is the shear stress of the velocity field, V the volume of the cell and Cs the Smagorinsky constant, which set to Cs ¼0.1. The factor fm takes the van Driest damping at the walls into account: e=AÞ3 , fm ¼ 1exp½ðy
ð2Þ
e ¼ rðywall u e : local veloe Þ=m (u with the normalized wall distance y city) and the damping factor A ¼25. lmax presents a cutoff in the case of volumes with large aspect ratios near the wall, i.e. with ywall 5 V 1=3 : lmax ¼ kywall . k was chosen to be 0.4. For all 3D calculations we used a hexahedral mesh with 931,000 volumes. A finer mesh was used at the boundaries up to a resolution of 0.12 mm at the first volume. Time steps have been chosen such that the root-mean-square Courant numbers were smaller than 1.
3. Characteristic numbers for convection The strength of the buoyancy convection is described by the Grashof number Gr ¼ g ar2 m2 l30 DT, where we chose the crucible
W. Miller et al. / Journal of Crystal Growth 318 (2011) 244–248
radius Rc as the characteristic length l0 and the maximum temperature in the melt (39 K) as DT. The other parametersp are ffiffiffiffiffiffigiven in Table 1. For analogy reasons one can identify Re ¼ Gr, where Re ¼ l0 u0 rm1 is the Reynolds number with a characteristic velocity u0. Thus, we obtain Rebuoy ¼ 1:8 104 :
ð3Þ
The Reynolds numbers for crucible and crystal rotation are defined by Rec ¼ 2pOc R2c rm1 and Rec ¼ 2pOX R2c rm1 , respectively. The rotation rates were OX ¼ 8:5 rpm and Oc ¼ 5 rpm for the crystal and crucible, respectively. With these values we obtain ReX ¼ 1:7 104 ,
Rec ¼ 9:8 103 :
ð4Þ
In the same way as for the buoyancy forces we can define a Reynolds number for the Lorentz forces as ReL ¼
maxðFL ÞR3c
rm2
:
ð5Þ
The maximum Lorentz force density max (FL) is about 2000 N m 3 (see below) and thus we get ReL ¼ 6 104 ,
ð6Þ
Table 1 Physical parameter used for the computations. Melting point temperature Density Specific heat Thermal conductivity Dynamic viscosity Thermal expansion factor Gravitational acceleration
Tm
1685 K 2420 kg m 3 1000 J kg 1 K 1 60 W m 1 K 1 7 10 4 Pa s 1.41 10 4 K 1 9.81 m s 2
r cP
l
m a g
245
which is significantly larger than any of the other Reynolds numbers. From u0 ¼ Rem=ðRc rÞ we can estimate the velocity range in a vertical plane of the melt, yielding u0 ¼ 230 mm/s. Indeed, we observed maximum velocities in this range for the run with TMF.
4. Results 4.1. Temperature field The experimental set-up with a heater magnet module (HMM) around the crucible consisting of three side coils and two bottom spirals is described in detail in [3]. For the numerical simulation we consider a growth situation, where the height of the melt was 41 mm (see Fig. 2). First of all, we computed the thermal field in the entire furnace using the axisymmetric code CrysMAS [4]. Applying the heater power as in the experiment we observed a very small radial temperature gradient. We like to note that this calculation did not include melt convection. The low temperature difference between crystal and crucible of about 9 K and in the entire melt of about 39 K is a result of the high thermal isolation. For the axisymmetric calculation we chose the maximal diameter of the crystal (see Fig. 1), i.e. 94 mm. The temperature data of CrysMAS are used as a boundary condition for the free melt surface in the 3D calculation. Because the temperature is almost the melting point temperature around the crystal it is not a large error when using the data of the axisymmetric simulation in 3D though the crystal is now not axisymmetric any longer. The temperature profile along the free interface of the melt can be seen in detail in Fig. 3. The crystal/ melt interface was computed to be slightly convex near to a flat interface. In the 3D calculations we use a completely flat crystal/ melt interface and neglect the wetting of crystal and crucible for simplifying the geometry and the grid generation. The crystal shape is taken as measured after the experiment (see Fig. 1): 94 mm diagonal diameter, 58 mm facet length. 4.2. Velocity field without TMF
94 mm
38 mm
152 mm
20 mm
crystal melt free surface
Fig. 1. Geometry of crystal and top of crucible as used for the computations. The data correspond to those of the crystal grown.
We start with 1/10 of the real Reynolds number and then increase the Reynolds numbers twice until the real values. For the runs with the real Reynolds numbers we use a time step of 0.01 s for the computation with the full Reynolds numbers. We observe the typical convection roll going up at the crucible’s side walls. In Fig. 7 (right) the averaged velocity field in a vertical cut is shown. Such a velocity field causes a heat transport from the sides of the crucible towards the crystal. From Fig. 6 (middle) it can be clearly seen that the radial temperature gradient near the free surface is increased. Because we fixed the temperature at the free surface this effect cannot be seen there. Despite this restriction, it is clear that maintaining a stable low radial thermal gradient would be not possible by crucible rotation. The other effect of the heat transport by the convection roll is the low temperature at the bottom of the crucible causing the danger of freezing (see Fig. 6).
+
+
+
+
39 mm
47 mm
74 mm
z = 30 mm
20 mm
z = 41 mm
z = 0 mm Fig. 2. Side view of geometry with locations of monitor points. The monitor point 2 is at the same distance from the axis as the {1 1 0} facet of the crystal, and point 3 is at the same distance as the diagonal section of the crystal. The monitor points are in a plane normal to the {1 1 0} facets.
-190
Fr / Nm−3
-1930
-1550
z
-1160
W. Miller et al. / Journal of Crystal Growth 318 (2011) 244–248
-770
246
r Fig. 3. Temperatures at crystal/melt interface and free melt interface as computed by CrysMAS. In the vicinity of the crystal the temperature is almost the melting point temperature, i.e. the radial gradient is nearly zero.
Fig. 5. Vectors of Lorentz force distribution and radial component Fr of the Lorentz force density.
1724
1715
1705
1695
1685
T [K]
Fig. 6. Comparison of the time-averaged (7 s) temperature field in a vertical plane for heat conduction only (top), buoyancy flow and rotation (middle), and buoyancy flow, rotation and TMF (bottom). Vertical plane is the same as in Fig. 7. Fig. 4. Instantaneous velocities on the free melt surface for the case without TMF. Thick arrows indicate the direction of crucible and crystal rotation.
At the free melt surface we can observe a recirculation of the flow behind the rotating edge of the crystal (see a snapshot in Fig. 4). The radially inwards and outwards directed flow is changing with time whereas the recirculation remains. The eddy viscosity is with a maximum value of msgs ¼ 1 104 Pa s significantly smaller than those of the melt, stating that the effects of numerically unresolved turbulent mixing is small. 4.3. Velocity field with TMF The Lorentz forces have been computed with CrysMAS and transferred to Ansys-cfx. The Lorentz forces are quite large with a maximum of about 2000 N m 3. Most important for the success of growing a Si crystal with {1 1 0} facets are the large radial components Fr in the lower side part of the crucible and the vertical components
at the bottom Fz directing upwards, i.e. opposite to the gravity (see Fig. 5). Such a Lorentz force distribution causes a convection roll at the crucible moving downwards at the side walls and upwards in the region of the crystal. The inversion of the convection roll by a TMF was demonstrated by Klein et al. in numerical simulations [8]. Here, we use much higher Lorentz forces yielding a Reynolds number more than four times larger than that of the buoyancy force and than that of the crucible rotation (see Section 3.2). Therefore, it can be expected that the flow induced by TMF is dominant. Indeed, the convection roll is completely inverted (see left side of Fig. 7) with a maximum time-averaged velocity of 115 mm/s compared to 15 mm/s for the case without TMF. Due to the higher velocities we had to use smaller time step (0.005 s) than in the case without TMF. Because the azimuthal velocities are now small compared to the vertical and horizontal velocities we did not observe a specific flow pattern due to the nearly quadratic shape of the crystal as in the case without TMF. The eddy viscosity is now up to msgs ¼ 4 103 Pa s in
W. Miller et al. / Journal of Crystal Growth 318 (2011) 244–248
247
Fig. 7. Time-averaged velocities (7 s) in a vertical plane with TMF (left half) and without TMF (right half). The lengths of the arrows correspond to the velocity but the arrow length is three times larger in the right half than in the left one. Maximum time-averaged velocity with TMF was 115 mm/s and without TMF 15 mm/s.
TMF 1690
TMF
1690
1688 1688
T [K]
1686 0
1
2
3
4
5
6
0
7
1
2
3
4
2
3
4
5
6
7
1706 1692 TMF
1702 1698
1690
1694 TMF 1690
1688 0
1
2
3
4
5
6
7 t [s]
0
1
5
6
7
Fig. 8. Temperature oscillations at the monitor points 1, 2, 3 and 4 as defined in Fig. 2. Case without TMF is red and case with TMF is blue as marked. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
certain regions, which is about six times the viscosity of the melt. A detailed analysis of the oscillating flow field and the distribution of the eddy viscosity will be published elsewhere. 4.4. Comparison of case without and with TMF We will compare the results without and with TMF in more detail. Firstly, we look on the time-averaged velocities in Fig. 7. In the case with TMF (left) we have a large convection roll fitting into the melt region below the free surface with a maximum velocity (time-averaged) of 115 mm/s. The flow beneath crystal is in the same direction as caused by crystal rotation and no shear flow occurs as in the case without TMF (right). In total, there is only a weak convection below the crystal. Without TMF we have the typical convection roll going upwards at the rim of the melt plus a counter vortex at the top outer part. As a result of the completely different velocity fields we obtain different temperature fields (see Fig. 6). The thermal gradient radial direction is much higher without TMF than with TMF, whereas the vertical gradient in the vicinity of the crystal/melt interface is higher in the case with TMF. The latter will stabilize a
flat crystal/melt interface. There is also no danger of freezing at the bottom in the case with TMF. The numerical calculations explain the experimental experience that the growth of a facetted crystal with a flat interface is possible using such a TMF. The time dependence at the four monitor points defined in Fig. 2 is given in Fig. 8 for both the case without TMF and that with TMF (marked). Oscillations with frequencies higher than 10 Hz have been suppressed by time-averaging. The following can be noticed:
The temperatures in the case with TMF are nearly the same at all monitor points.
The frequency is higher in the case with TMF, most pronounced at monitor point 1.
The oscillation amplitude is smaller in the case with TMF. 5. Conclusions We have performed 3D simulations of the melt flow with and without TMF. The TMF as used in the experiments causes a strong vortex with a heat transport from the crucible’s bottom to the crystal. Such a heat transport stabilizes the low temperature gradient along
248
W. Miller et al. / Journal of Crystal Growth 318 (2011) 244–248
the free surface and the stable growth of the silicon crystal with {1 1 0} facets becomes possible. The TMF induced vortex is in opposite direction to the common vortex caused by buoyancy convection and enhanced by crucible rotation. The results of the numerical calculations show that this rotation-driven vortex leads to a significant increase of the radial gradient near the free surface, which would inhibit the formation of {1 1 0} facets because of the decreased undercooled regions along the {1 1 0} places (see [3]). References [1] A. Muizˇnieks, A. Rudevics, K. Lacis, H. Riemann, A. Ludge, F.W. Schulze, B. Nacke, Square-shaped silicon crystal rod growth by the fz method with special 3D shaped HF inductors, Magnetohydrodynamics 43 (2) (2007) 269–282. [2] P. Rudolph, Travelling magnetic fields applied to bulk crystal growth from the melt: the step from basic research to industrial scale, J. Crystal Growth 310 (2008) 1298–1306.
[3] P. Rudolph, M. Czupalla, B. Lux, F. Kirscht, C. Frank-Rotsch, W. Miller, M. Albrecht, The use of heater-magnet module for Czochralski growth of PV silicon crystals with quadratic cross section, J. Crystal Growth, this issue, doi:10.1016/j.jcrysgro.2010.10.070. [4] J. Fainberg, D. Vizman, J. Friedrich, G. Mueller, A new hybrid method for the global modeling of convection in CZ crystal growth configurations, J. Crystal Growth 303 (2007) 124–134. [5] W. Miller, Numerical simulations of bulk crystal growth on different scales: silicon and GeSi, Phys. Status Solidi B 257 (4) (2010) 855–869. [6] A. Muiznieks, A. Krauze, B. Nacke, Convective phenomena in large melts including magnetic fields, J. Crystal Growth 303 (2007) 211–220. [7] I.Y. Evstratov, V.V. Kalaev, A.I. Zhmakin, Y.N. Makarov, A.G. Abramov, N.G. Ivanov, A.B. Korsakov, E.M. Smirnov, E. Dornberger, J. Virbulis, E. Tomzig, W. von Ammon, Numerical study of 3D unsteady melt convection during industrial-scale CZ Si-crystal growth, J. Crystal Growth 237–239 (3) (2002) 1757–1761. [8] O. Klein, C. Lechner, P.-E. Druet, P. Philip, J. Sprekels, C. Frank-Rotsch, F.-M. Kießling, W. Miller, U. Rehse, P. Rudolph, Numerical simulation of Czochralski crystal growth under the influence of a traveling magnetic field generated by an internal heater-magnet module (HMM), J. Crystal Growth 310 (2008) 1523–1532.