ARTICLE IN PRESS
Journal of Crystal Growth 266 (2004) 140–144
A three-dimensional numerical simulation study of the Marangoni convection occurring in the crystal growth of SixGe1x by the float-zone technique in zero gravity Hisashi Minakuchia, Yasunori Okanoa,*, Sadik Dostb a
Department of Materials Science and Chemical Engineering, Shizuoka University, Johoku 3-5-1, Hamamatsu 432-8561, Japan b Crystal Growth Laboratory, University of Victoria, Victoria, BC, Canada V8W 3P6
Abstract A numerical simulation study was carried out to examine the effect of solutal Marangoni convection during the SixGe1x crystal growth by the floating-zone (FZ) technique under zero gravity. In this study, a three-dimensional fullzone configuration for the FZ was considered. Computational results show that the contribution of the solutal Marangoni convection to the flow structure was significant although the strength of the solutal Marangoni convection is much weaker than that of the thermal Marangoni convection. Application of a crystal and feed rotation with an optimum rate induces a good mixing in the melt and leads to crystals with a uniform concentration distribution along the growth interface. r 2004 Elsevier B.V. All rights reserved. PACS: 02.60.Cb; 02.70.Bf; 47.20.Dr; 47.27.Te Keywords: A1. Computer simulation; A1. Fluid flow; A2. Floating zone technique; B1. Germanium silicon alloys
1. Introduction The alloys of silicon (Si) and germanium (Ge) (i.e., SixGe1x) are promising semiconductors materials for future novel opto-electronic devices that will be used in various fields of the information technology as optical switches, photo-detectors, solar cells, and power generators due to their low cost, light weight, stable performance and low energy requirement. In all these applications, the compositional uniformity of grown crystals used in the form of wafers and substrates is of utmost *Corresponding author. Tel./fax: +81-53-478-1169. E-mail address:
[email protected] (Y. Okano).
importance. The growth of SixGe1x crystals with uniform composition (x) is not easy. This is mainly due to the large gap between the liquidus and solidus curves in the phase diagram. The gravitational field of Earth makes the situation worse because of the significant density difference between Ge and Si; the density of Ge is much larger than that of Si, ðrGe :5:51 g=cm3 > rSi : 2:52 g=cm3 Þ leading to further segregation in the melt. The floating zone (FZ) method in a microgravity environment may provide a great platform for the growth of high-quality crystals. This is mainly due to two reasons, namely, (i) in this technique there is no possibility of crucible contamination, and (ii) there is less gravitational segregation due to a weaker
0022-0248/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2004.02.038
ARTICLE IN PRESS H. Minakuchi et al. / Journal of Crystal Growth 266 (2004) 140–144
gravitational force in such an environment. However, even in a microgravity environment, the unsteady and three-dimensional Marangoni convection in the melt incorporates nonuniformity into the grown crystals. In this respect, a good understanding for the relative contribution of Marangoni convection is essential. From this point of view, a large number of three dimensional analyses have been carried out for the thermal Marangoni convection in FZ systems (see, for instance,Refs. [1–6]). In Ref. [6], the authors considered an actual FZ system and performed a three-dimensional numerical simulation for the thermal Marangoni convection in the FZ. It was shown that the flow and concentration fields become quickly three dimensional, in spite of assumed axisymmetric thermal boundary conditions, although the thermal field remained almost twodimensional during the computational time frame. However, in the growth of alloy SixGe1x crystals, not only the thermal Marangoni convection but also the solutal Marangoni convection plays a significant role [7], and leads to further segregation in the melt. The present study aims at the understanding of both the thermal and solutal Marangoni flows during the crystal growth of SixGe1x by the FZ method. In order to identify clearly the role of Marangoni convection, the contribution of the buoyancy convection was removed by considering zero gravity in the model. The effect of crystal rotation was also included in the study to determine its effect on concentration uniformity along the growth interface.
2. Numerical method Fig. 1 shows the schematic configuration of the model in which a liquid bridge between two rods of the same diameter is heated from the outside perimeter. The ambient nondimensionl temperature profile, YaðzÞ; which is defined as ðTa T0 Þ=ðTH Tc Þ; is represented by the following equation [8]: ( ) 1 2 ; YaðzÞ ¼ YC þ ðYH YC Þ exp 25 Z 2 ð1Þ
141
Si0.05Ge0.95(poly) cold disc ΘC ΘH ΘC 0
Z=1
L
R=1
ΘH > ΘC SiGe(single) cold disc Fig. 1. Schematic configuration of the model.
where YC and YH are nondimensionl minimum and maximum temperatures, respectively. In the model the following assumptions were made. The liquid zone is taken as an incompressible, Newtonian fluid. In an actual FZ growth system, the feed/melt (melt) and the seed/melt (growth) interfaces will not be flat; will follow the isotherms according to the phase diagram. For computational convenience and simplicity, however, the deformation of the melt and the growth interfaces (solid/liquid interfaces) on the flow structures are neglected, and then the interfaces considered as flat surfaces in the model. In order to examine the individual role of the Marangoni convection in the liquid zone, the gravitational constant was taken zero in the model. This assumption eliminated the contribution of the buoyancy convection in the melt, and also made possible to consider the liquid/gas interface as a flat cylindrical surface. This assumption may not be a reality in any microgravity environment such as the Space Station or Space Shuttle due to the presence of a time-depended residual gravity field, but it is a good approximation for examining the role of the Marangoni convection in the FZ melt, and also for determining the relative contribution of crystal rotation on the transport structures of the liquid zone. Finally, the heat transfer between the liquid zone and the ambient is considered to be in all forms, i.e., conductive, convective and radiative. The composition of the polycrystalline material is selected as Si0.05Ge0.95. Under these assumptions, the continuity, momentum, energy and mass transport equations,
ARTICLE IN PRESS H. Minakuchi et al. / Journal of Crystal Growth 266 (2004) 140–144
142
and the associated boundary conditions are presented in the three-dimensional cylindrical coordinate system (these equations are not given here for brevity; details can be found in Ref. [6]). Boundary conditions are given as follows: Along the solid/liquid interfaces (Z ¼ 0 and 1): U ¼ 0; Y ¼ 0:5; V ¼ Vp ; W ¼ Reo R ðat Z ¼ 0Þ and W ¼ Reo R ðat Z ¼ 1Þ; ð2Þ ðC Cf ÞVp ¼ Sc
qC ðZ ¼ 1Þ; qZ
Cð1 k0 ÞVp ¼ Sc
ð3Þ
qC ðZ ¼ 0Þ: qZ
ð4Þ
Along the free surface (R ¼ 1): qY ¼ ðRa0 G þ BiÞðY YaðZÞÞ; qR W qY qC 2 q þ MaC ; ¼ MaT R qR R qy qy qV qY qC ¼ MaT þ MaC : qR qZ qZ
ð5Þ ð6Þ
ð7Þ
In these equations, U; V ; W ; C; Cf ; Vp ; Reo ; Bi; Ra0 and G are used respectively to denote nondimensional radial velocity (Lu=n), nondimensional axial velocity (Lv=n), nondimensional azimuthal velocity (Lw=n), molar fraction of silicon, molar fraction of silicon in the polycrystal (Cf ¼ 0:05), nondimensional growth velocity, rotational Reynolds number, Biot number, Ra0 ¼ Less =l and G ¼ ½T 2 þ Ta2 ðzÞ½T þ TaðzÞ; respectively. Eqs. (6) and (7) are the tangential force balance at the free surface which describes the contribution of both thermal Marangoni convective and solutal Marangoni convective flows. The associated dimensionless thermal Marangoni and solutal Marangoni numbers are defined as the following: MaT ¼
DTmax L qs ; mn qT
MaC ¼
The governing equations and the boundary conditions were discretized by the finite difference method on a staggered grid system, and solved by the HSMAC method [6]. After examining mesh dependency [6], the numerical domain was divided into NR Ny NZ ¼ 31 31 31 mesh with uniform spacing. The third-order upwind difference schemes and the fourth-order central difference schemes were applied, respectively, to the inertia and viscous terms in the govern ing equations. Physical properties of silicon– germanium melt were used (Pr ¼ 7:71 103 ; Sc ¼ 22:5 and k0 ¼ 5; where Pr; Sc and k0 stand for the Prandtl number, Schmidt number and segregation coefficient, respectively) [7]. The length of free surface, aspect ratio and growth velocity were set to L ¼ 1:0 102 m, Asp ¼ 1 and vp ¼ 1:39 107 m/s, respectively. The results of the conductive heat transfer analysis were used as the initial conditions of computations.
3. Results and discussion Fig. 2 represents computational results in the absence of thermal Marangoni convection. In this figure, the left and right halves describe the concentration and flow fields in the vertical plane, respectively. As can be seen, the flow in the liquid zone is in the direction from the lower disk to upper disk along the free surface. This flow pattern is similar to that of the thermal Marangoni convection in a Half-zone configuration. Fig. 3 shows the time dependency of maximum flow strength ðjV jMax Þ: The strength of the induced flow
DCmax L qs ; mn qC
where DTmax is the maximum temperature difference, DCmax the maximum silicon concentration difference, m the viscosity, n the kinematic viscosity and s the surface tension.
Fig. 2. Concentration (left) and flow (right) fields at MaT ¼ 0 and Mac ¼ 1 105 (t ¼ 150 s, DC ¼ 0:5 103 ; DC is the difference of silicon concentration).
ARTICLE IN PRESS H. Minakuchi et al. / Journal of Crystal Growth 266 (2004) 140–144
143
50
|V|Max[-]
40 30 20 10 0
0
50
100
150
200
Time [sec] Fig. 3. The time dependency of maximum flow strength ðjV jMax Þ when MaT ¼ 0 and Mac ¼ 1 105 :
Fig. 6. The silicon concentration field along the free surface when MaT ¼ 5 103 and Mac ¼ 1 105 at t ¼ 140 s (DC ¼ 0:5 103 ).
|V|Max [-]
600 550 500
Ma ΜΑ TΤ C =0 3 5 Ma T =5×10 ΜΑ ΧΤ ,Ma C =1×10 =5×10 3 ,Ma
450 400
0
100
200
300
400
Time [sec] Fig. 4. The time dependency of maximum flow strength ðjV jMax Þ when the thermal Marangoni convection is present.
Fig. 5. Time dependency of the iso-velocity when MaT ¼ 5 103 ; Mac ¼ 0 (a) and MaT ¼ 5 103 ; (b) Mac ¼ 1 105 :
(solutal) velocity is much smaller than that by the thermal Marangoni convection in spite of the fact that the solutal Marangoni number was assumed to be larger than the thermal Marangoni number. Figs. 4 and 5 show, respectively, the time dependency of maximum flow strength ðjV jMax Þ and the iso-velocity face when the thermal Marangoni convection is also present. In early stages of the simulation, the flow field is still two
dimensional, and two vortexes of the same size develop in the melt. As time proceeds, however, the flow field becomes three-dimensional. We note that when both the solutal and thermal Marangoni convection exist, the flow becomes three dimensional earlier than that of the thermal Marangoni convection alone. As can been seen from Fig. 5, the solutal Marangoni convection significantly affects the flow structure although the contribution of the solutal Marangoni convection was very small as discussed in Figs. 2 and 3. The structures of the solutal Marangoni convection change with m ¼ 4 to that with m ¼ 3: Fig. 6 shows the concentration distribution in the y-direction, which leads to the solutal Marangoni convection in the y-direction. From the above results, we can conclude that the solutal Marangoni convection in the y-direction significantly contributes to the flow structures and the concentration distribution in the melt although the strength of the induced velocity due to the solutal Marangoni convection is weak. In Ref. [6], the authors showed that the rotation of the crystal and the feed was beneficial in growing axisymetrically uniform crystals. However, only one rotation, 5 rpm, was considered in Ref. [6]. In the present study, in order to have a better idea for the effect of crystal rotation, two additional rotation rates, 10 and 15 rpm, were applied. The effect of crystal rotation on the concentration distribution along the growth interface was discussed. Based on the numerical results at 500 s at MaT ¼ 5 103 and Mac ¼ 1 105 ; it
ARTICLE IN PRESS H. Minakuchi et al. / Journal of Crystal Growth 266 (2004) 140–144
144
15 rpm (Reo ¼ 1:16 103 ). This is because the rotation induced flow at this rate was dominant. One can therefore conclude that a rotation rate of 15 rpm is very beneficial since it leads to a better mixing in the melt.
Fig. 7. The distribution of the silicon concentration along the growth interface when MaT ¼ 5 103 and Mac ¼ 1 105 at t ¼ 530 s (DC ¼ 0:5 103 ).
0rpm 5rpm
10rpm 15rpm
0.004
dCMax [-]
0.003 0.002 0.001 0.000 450
500
550
600
Time [sec] Fig. 8. The time dependency of maximum silicon concentration difference in the growth interface when MaT ¼ 5 103 and Mac ¼ 1 105 :
was assumed that both the upper and lower rods rotate at the same speed but in the opposite directions. Computations were repeated in such a case. Fig. 7 shows the effect of rotation on the silicon concentration distribution along the growth interface at t ¼ 530 s. The time dependency of the maximum silicon concentration difference (dCMAX ) in the growth interface is presented in Fig. 8. As can been seen from these figure, the application of a crystal rotation is beneficial in growing axisymetrically uniform crystals. However, low rotation rates, for instance 5 rpm, induce a concentration distribution of ‘‘W ’’ shape as pointed out in Ref. [6]. As shown in Fig. 8, the values of dCMAX at 5 rpm (Reo ¼ 3:88 102 ) and 10 rpm (Reo ¼ 7:75 102 ) are almost the same. On the other hand, dCMAX becomes very small at
4. Conclusion The effects of solutal Marangoni convection and crystal rotation were numerically studied. A threedimensional floating-full-zone liquid bridge for the growth of Si/Ge under zero gravity was considered. The study led to the following results: 1. Strength of the axial flow along the free surface due to the solutal Marangoni convection was much weaker than that of the thermal Marangoni convection although the assumed solutal Marangoni number was larger than the thermal Marangoni number within the limit of the present computations. 2. The account for the solutal Marangoni convection changes the evolution of the flow in time. 3. Rotation of the crystal and the feed was very beneficial in growing axisymmetrically uniform crystals. In order to obtain a more uniform concentration distribution along the growth interface, an optimum rotation rate that induces a good mixing in the melt must be used.
References [1] H.C. Kuhlmann, Thermocapillary Convection in Models of Crystal Growth, Springer Tracts in Modern Physics, Vol. 152, Springer, Berlin, Heidelberg, 1999. [2] Th. Kaiser, K.W. Benz, J. Crystal Growth 183 (1998) 564. . M. Lichtensteiger, Th. Kaiser, K.W. [3] P. Dold, A. Croll, Benz, J. Crystal Growth 231 (2001) 95. [4] C.W. Lan, J.H. Chian, J. Crystal Growth 230 (2001) 172. [5] N. Imaishi, S. Yasuhiro, Y. Akiyama, S. Yoda, J. Crystal Growth 230 (2001) 164. [6] H. Minakuchi, Y. Okano, S. Dost, Int. J. Mater. Product. Tech., accepted. . K.W. Benz, [7] T.A. Campbell, M. Shweizer, P. Dold, A. Croll, J. Crystal Growth 226 (2001) 231. [8] Y. Okano, N. Audet, S. Dost, S. Kunikata, Chem. Eng. J. 84 (2001) 315.