Three-dimensional instabilities in Czochralski process of crystal growth from silicon melt

ARTICLE IN PRESS

Journal of Crystal Growth 305 (2007) 185–191 www.elsevier.com/locate/jcrysgro

Three-dimensional instabilities in Czochralski process of crystal growth from silicon melt Y. Rosenstein, P.Z. Bar-Yoseph Computational Mechanics Laboratory, Faculty of Mechanical Engineering, Technion-Israel Institute of Technology, Israel Received 17 May 2006; received in revised form 7 January 2007; accepted 29 March 2007 Communicated by J.J. Derby Available online 20 April 2007

Abstract This paper deals with axisymmetry breaking instabilities in Czochralski process of crystal growth from silicon melt (Prandtl number Pr ¼ 0:01). Numerical 3D linear stability analysis was carried out on the axisymmetric bulk flow model. Stability diagram of critical Grashof numbers Grc dependent on aspect ratio a( ¼ height/radius) in the range 0:4pap1:0 was computed. Computations were carried out using the spectral element method in the meridional plane with Fourier decomposition in the azimuthal direction. It was found that convective instability sets in through an Hopf bifurcation displaying oscillations in time. Sensitivity of mode transitions was observed at parameter range of a40:65 and in some regions modes were observed approaching each other closely. Dangerous modes in the 0:65oap0:8 region were 2, 3, 4. For 0:4pap0:85 dispersion relation analysis reveals convective instability effects while for larger a rotational effects appear. Analysis for the case of no heat convection ðPr ¼ 0Þ was carried out showing different behavior for 0:4pap0:85, thus, validating the analysis of Pr ¼ 0:01. r 2007 Elsevier B.V. All rights reserved. PACS: 81.10; 47.20 Keywords: A1. Computer simulation; A1. Convection; A1. Crystal growth; A2. Czochralski method

1. Introduction Czochralski based crystal growth processes may display transitions from steady axisymmetric flow into asymmetric time-dependent flows. This in turn may cause inhomogeneities in the grown crystal [1,2]. The dynamics of the flow are complex and require the solution of the 3D time-dependent flow equations coupled with the heat equation. Many published works address the direct simulation of timedependent 3D problem (e.g. Ref. [3]); however, this is CPU-time consuming and requires the simultaneous solution of millions of equations depending on an infinite manifold of initial conditions and hence is not general. An alternative approach employing the methods of hydrodynamic stability analysis [4,5], offers considerable reduction in computer resource usage and is general to all Corresponding author. Tel.: +972 4 8293476.

E-mail address: [email protected] (Y. Rosenstein). 0022-0248/$ - see front matter r 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2007.03.057

infinitesimal perturbations. In this approach, the stability of axisymmetric steady flow with respect to 3D timedependent perturbations is analyzed and critical Grashof number and frequency are computed. This approach is advantageous in crystal growth processes, since it enhances the understanding of the transition from steady axisymmetric regime to 3D time-dependent flows. The 3D timedependent problem is converted into a sequence of axisymmetric problems in Fourier space, simplifying the numerical treatment of the problem and saving computer resources. Similar works have shown relevance to CZ 3D flows for specific Prandtl number Pr ¼ 1:4 (see Ref. [6]). This work approaches the problem using bulk flow modeling [7] based on the international test [8]. The spectral elements method pioneered by Patera [9] is used to discretize the steady axisymmetric Navier–Stokes equations coupled with the equation of energy through the Boussinesq approximation. Pressure is eliminated using the consistent penalty method [10]. The equations are then

ARTICLE IN PRESS 186

Y. Rosenstein, P.Z. Bar-Yoseph / Journal of Crystal Growth 305 (2007) 185–191

assembled and solved using preconditioned GMRES method [11]. Three-dimensional time-dependent perturbations are superimposed on the steady solution using Fourier decomposition of the azimuthal direction. The resulting linear eigenvalue problem is then solved using subspace iterations [12,13]. This paper is organized as follows: In Sections 2 and 3 we briefly describe the mathematical model and solution technique employed in this work. In Section 4 results are displayed and described for silicon melt ðPr ¼ 0:01Þ. Section 5 concludes with a brief summary.

The prescribed boundary conditions are (Fig. 1): qT ¼0 qz on z ¼ 0. T ¼1

u ¼ 0;

We consider a co-axial cylinder-disk configuration [7,8], where the disk represents the crystal, the disk is free to rotate and the crucible is held fixed (Fig. 1). We employ the standard Navier–Stokes equations coupled with the energy equation through the Bossinesq approximation. We denote by Rc , Rx , T c , T x , Ox the crucible and crystal radii, temperatures and crystal’s angular velocity, respectively. By r, n, k, cp , g we denote the density, kinematic viscosity, thermal conductivity, heat capacity and coefficient of thermal expansion of the melt, respectively. Length, velocity and temperature are normalized by Rc , n=Rc and ðT c  T x Þ, respectively. Thus, the non-dimensional equations describing the flow are

(5)

at r ¼ 1. ur ¼ uz ¼ 0;

uy ¼ rRex ;

T ¼0

(6)

rb 1b

(7)

at 0prpb; z ¼ a. qur quy ¼ ¼ uz ¼ 0; qz qz

2. Mathematical formulation of the problem

(4)

u ¼ 0;

T¼

at borp1; z ¼ a. Pole conditions are defined at r ¼ 0 (see Refs. [14,15]): ur ¼ uy ¼

qvz qT ¼ 0; ¼ qr qr

qur quy ¼ ¼ uz ¼ T ¼ 0; qr qr ur ¼ uy ¼ uz ¼ T ¼ 0;

m ¼ 0,

(8)

jmj ¼ 1,

(9)

jmj41,

(10)

m is the Fourier wave number defined in Eq. (17). Periodicity is assumed in the azimuthal direction. The dimensionless parameters appearing in Eqs. (1)–(10) are defined as follows: a¼

H Aspect ratio (height/crucible radius), Rc

(11)

(2)

b¼

Rx Ratio of crystal to crucible radii, Rc

(12)

(3)

Rex ¼

R2c Ox Crystal Reynolds number, n

(13)

qu þ ðu  rÞu ¼ rp þ r2 u þ Gr Tez , qt

(1)

r  u ¼ 0, qT 1 2 þ ðu  rÞT ¼ r T, qt Pr

where u ¼ ður ; uy ; uz Þ denote the velocity components in cylindrical coordinates.

Gr ¼

ggðT c  T x ÞR3c Grashof number, n2

(14)

Pr ¼

nrcp Prandtl number. k

(15)

3. Method of solution Two methods of analysis exist to solve the problem, direct 3D numerical solution or 3D linear stability analysis on the basic (axisymmetric) solution. The direct method requires discretization in 3D space, time marching and infinite variety of initial conditions, this renders the problem very difficult to solve and the method is not general. Instead we choose the method of 3D linear stability analysis on the axisymmetric base flow, which is general to all infinitesimal perturbations. Denoting by x ¼ ður ; uy ; uz ; TÞ Fig. 1. Czochralski process of crystal growth [8]—sketch of the problem.

(16)

the vector of velocities and temperature, x is written as superposition of axisymmetric steady flow and 3D

ARTICLE IN PRESS Y. Rosenstein, P.Z. Bar-Yoseph / Journal of Crystal Growth 305 (2007) 185–191

time-dependent perturbations: x ¼ Xðr; zÞ þ

1 X

xm ðr; zÞelm tþimy .

187

found from (17)

Grc ¼ min Grc ðmÞ. m

(22)

m¼1

Substituting Eq. (17) into Eqs. (1)–(10) the problem is now decomposed into a combination of nonlinear steady axisymmetric problem and 3D linear eigenvalue problems. We employ the spectral element method to discretize the axisymmetric components of the equations, using penalty method to eliminate pressure with penalty parameter  ¼ 107 . Thus the nonlinear (axisymmetric) part of the problem can be written as ðA þ CðuÞÞu þ Bp  Gr MTez ¼ 0,

(18)

ðA þ CðuÞÞT ¼ 0,

(19)

BT u ¼ 0,

(20)

where velocity space is taken to be PN (Polynomials up to order N) and pressure space is taken as PN2 in order to satisfy the LBB condition [16]. Boldfaced operators A and C are the vector diffusion and convection operators, respectively. A and C represent their scalar counterparts, respectively. M is the discrete lumped mass operator, B is the discrete gradient operator. The system is then linearized using Newton’s method and arc-length continuation is employed to march on different solution branches. Solution of the linearized system is obtained using preconditioned GMRES method. Thus the basic solution is obtained. Substituting the 3D perturbations into Eqs. (1)–(10) and neglecting second order terms in perturbations we obtain the generalized eigenvalue problem: Gx ¼ lm Mx,

(21)

where the G operator represents the linearized perturbation equations and M is again the spectral mass matrix. This problem is solved using the method of subspace iterations. The eigenvalue lm is in general complex, lm ¼ sm þ iom with sm and om real. When sm 40 the flow is unstable. Stability analysis is performed in the following way: (i) given parameters ðGr; Pr; a; b; Rex Þ, loop over the next steps; (ii) given m, compute basic steady axisymmetric solution to Eqs. (18)–(20); (iii) compute eigenvalue problem (21) to obtain eigenvalue lm and eigenvector; (iv) increase Gr and repeat steps (i), (ii), (iii) until the real part of lm changes sign; (v) when the real part of lm changes sign Grc ðmÞ is obtained. Thus Grc ðmÞ, the critical Gr number for each Fourier mode m is found. The dominant critical values are

4. Results In all our computations the polynomial order of each element was maintained at 7  7. 4.1. Code validation The eigensolver was compared with Ref. [17], onset of convection in cylindrical cavity. A cylindrical cavity of unit aspect ratio is heated from below with anti-symmetric temperature boundary conditions at top and bottom. The steady solution is conduction. Dominant Fourier modes and Grc were compared. In Ref. [17] direct numerical simulation of the time-dependent Navier–Stokes equations using spatial discretization of five spectral elements of 7  7  9 nodal points each for the x; y; z directions, respectively, was carried out. We compared the results of our linear stability analysis for modes 0, 1, 2 with those of Ref. [17]. The mesh was taken at 7  7 elements. Comparison of our results with those of Ref. [17] show very good agreement for modes 0; 1; 2 (Grc ¼ 36160; 38940; 41783, respectively) with largest relative error of 0:4% for the case m ¼ 1. 4.2. Czochralski process Simulations were carried out on the international test problem (see Fig. 1) in the range of parameters Pr ¼ 0:01, 0:4pap1:0. All simulations were carried out with the constant parameter values Rex ¼ 500 and b ¼ 0:4. For silicon melt n ¼ 3:1  107 m2 =s, r ¼ 2750 kg=m3 , Pr ¼ 0:01. Rex ¼ 500 represents rotation rate of 0:6 RPM for a 50 mm radius crucible, 20 mm radius disk. a ranges represent stages in the process. Eigensolver convergence was tested for the two most dangerous modes m ¼ 2; 3 at a ¼ 0:65 using 5  5, 7  7 and 9  9 elements. It was found that at Gr ¼ 1:3  106 both modes 2 and 3 were dominant. Based on the convergence tests, resolution of 7  7 elements was chosen. Typical steady axisymmetric picture for the range of a considered can be seen in Fig. 2 for a ¼ 0:6. Figs. 3 and 4 display the dependence on a of the first 5 modes for critical Grc ðmÞ and oc ðmÞ, respectively. The critical curves show that at ao0:65 mode 3 is dominant while many mode transitions occur at aX0:65 which are accompanied by either sharp changes in critical Gr or critical o. However, the changes in Grc and oc do not have to be coincident. Also modes approach each other very closely at 0:65pap0:8. There is evidently, sensitivity of dominant modes to aspect ratio. Modes 2, 3, 4 compete closely for dominance at a ¼ 0:65. At a ¼ 0:7 modes 2, 4 compete with

ARTICLE IN PRESS Y. Rosenstein, P.Z. Bar-Yoseph / Journal of Crystal Growth 305 (2007) 185–191

188

0.6 0.5 0.4

0.4

0.3 0.2

0.2 0.1 0.2

0.4

0.6

0

0.8

0

0.2

0.4

0.6

0.8

Fig. 2. Streamlines (left) and isotherms (right) for a ¼ 0:6, Gr ¼ 1:3  106 .

x 103

3

m=0 m=1 m=2 m=3 m=4

2.5 2

10

1.5

Grc

Grc (m)

x 105

15

m=1

1

m=2

5

m=3 m=4

0.5 0

0.4

0.5

IV

III

II

l

0.6

0.7 α

0.8

0.9

0

1

0.4

Fig. 3. Critical GrðmÞ numbers vs. a with wave number as parameter. 1400 III

II

I

IV

1200

0.5

0.6

0.7 α

0.8

0.9

1

Fig. 5. Stability diagram showing critical Gr numbers vs. a.

103

m=0 m=1 m=2 m=3 m=4

I

II

III

IV

800

ωc

ωc(m)

1000

600

102

400

m=1 m=2 m=3 m=4

200 0 0.4

0.5

0.6

0.7 α

0.8

0.9

1

101 0.4

0.5

0.6

α

100

0.8

Fig. 4. Critical frequencies am vs. a with wave number as parameter.

Fig. 6. Log–log plot of the critical frequencies as function of a.

each other closely. Sensitivity of modes to aspect ratio can also be seen in the range 0:8pap1:0 although not all modes approach each other closely. The resulting stability diagram which is found from Eq. (22) is displayed in Fig. 5.

To obtain some quantitative analysis let us define the dispersion relation as function of m; a; Pr (see for example Ref. [18]): oc ¼ f ðm; a; Pr; Rex Þ.

(23)

ARTICLE IN PRESS Y. Rosenstein, P.Z. Bar-Yoseph / Journal of Crystal Growth 305 (2007) 185–191

A log–log plot of the curve is depicted in Fig. 6, it is obtained by taking the frequency of the most dominant mode as function of a. It is clearly seen that at Sections 1 and 2 the curve is almost linear in logðaÞ with mean slope 4. In Section 3 the curve is almost constant. In Section 4 the curve is oscillatory. Thus, at Sections 1 and 2 one can deduce that oc  Oða4 Þ, at Section 3 o ¼ f ðmÞ with no dependence on a and at Section 4, the behavior is oscillatory. The behavior of oc in Sections 1–3 is typical of convective instability (see for example Ref. [4]), while the behavior of oc at Section 4 is not typical of convective instability. Based on the previous discussion we hypothesize that two different mechanisms dominate the instability depending on a. To establish that, further analysis was carried out for Pr ¼ 0 the results of which are shown in

6

m=0 m=4

5 4 Grc

Figs. 7 and 8. What is clearly seen is anomalous behavior of the Pr ¼ 0 stability curve. For ao0:85 axisymmetric mode ðm ¼ 0Þ is dominant with steady convection. For 0:85pap1 mode 4 is dominant with small frequencies. This indicates that in the absence of temperature convection the hydrodynamic instability mechanism is steady axisymmetric up to a ¼ 0:85, from a40:85 the dominant mode is 4 with weak oscillatory behavior. We can now deduce that at a40:85 rotational effects play a dominant part in the instability mechanism. From the foregoing discussion it is apparent that in the range of geometric and rotational parameters considered for Pr40 the dominating convection mechanism is mixed both hydrodynamic and thermal. At a40:85 rotational effects become dominant.

4.3. Case study: comparison of stability analysis with time-dependent 3D works

x 105

3 2 1 0 0.4

189

0.5

0.6

0.7 α

0.8

0.9

1

Fig. 7. Stability diagram showing critical Gr numbers for Pr ¼ 0.

In this section, we show that results of our stability analysis correspond to results obtained by 3D timedependent simulations. Comparison is made with the work of Tomonari et al. [19] who have computed the 3D timedependent flow of molten silicon (Pr ¼ 0:016 in their work) in shallow Czochralski configuration (a ¼ 0:5, b ¼ 0:25) at constant Ra ¼ 105 and 106 . It was found in that work that at Ra ¼ 105 the axisymmetric flow breaks into a mode 3 asymmetric flow consisting of six standing rolls forming a triangular shape. We have performed linear stability analysis at Ra ¼ 105 ðGr ¼ 6:25  106 Þ. It was found that mode 3 was dominant with oscillation frequency of o ¼ 0:21. This indicates that it is not standing but oscillating in very low frequency and therefore may appear standing. The basic temperature superposed with the mode 3 perturbation is plotted in Figs. 9 and 10 for two different heights. 1

0.03

0.8

m=0 m=4

0.6

0.025

0.4

ωc

0.02

0.2 0

0.015

−0.2

0.01

−0.4 −0.6

0.005

−0.8

0

0.4

0.5

0.6

0.7 α

0.8

Fig. 8. Critical frequencies for Pr ¼ 0.

0.9

1

−1

−1

−0.5

0

0.5

1

Fig. 9. Perturbed temperature distribution at z ¼ 0:9a for Gr ¼ 3:0  106 , Pr ¼ 0:016, a ¼ 0:5, b ¼ 0:25, Rex ¼ 0.

ARTICLE IN PRESS Y. Rosenstein, P.Z. Bar-Yoseph / Journal of Crystal Growth 305 (2007) 185–191

190 1

terms of the Fourier series composed from the perturbations (multiplied by a proper amplitude). Our analysis shows that a finite number of terms is required to reconstruct the flow. Our results reveal that was seen as standing waves in the numerical and experimental work of Tomonari et al. [19] is actually a very slowly rotating wave. Very slowly changing phenomena are hard to measure or compute in direct numerical simulation. Stability analysis is thus a complementary tool which enables the observation of these phenomena [21]. Numerical convergence was tested and shows convergence of critical numbers at relatively coarse mesh of 7  7 elements of 7  7 polynomial order each.

0.8 0.6 0.4 0.2 0 −0.2 −0.4 −0.6 −0.8 −1

−1

−0.5

0

0.5

Fig. 10. Perturbed temperature distribution at Gr ¼ 3:0  106 , Pr ¼ 0:016, a ¼ 0:5, b ¼ 0:25, Rex ¼ 0.

Acknowledgments

1

z ¼ 0:08a

for

Thus, we were able to reconstruct the 3D flow using the superposition of basic flow and Fourier dominant terms. An argument was made in Ref. [19] that the radiation condition on the free surface is the cause for the transition to asymmetry. However, our computations show that this transition is inherent also in the bulk flow model that we employed without radiation. This indicates that more insight into the flow structure can be gained by stability analysis. 5. Concluding remarks The present paper reports preliminary results of the study of three-dimensional instabilities of an axisymmetric flow model of Czochralski process. Further work would include the effects of Marangoni flow and would have to include a kinetic model of the mass transfer. Future work should include the effects of electromagnetic fields as well, since it has been shown [20] that these have stabilizing effect on the flow. It should be emphasized that this model was chosen for its simplicity and yet its spectral behavior is rich. Much insight is gained from the results obtained from this model on the instabilities. It is shown that the destabilizing mechanism is mixed thermal and rotational. Modes turned out to be sensitive to changes in geometry (aspect ratio). For aspect ratio range between 0.65 and 0.8 competing modes were 2, 3, 4. From the dispersion relation analysis it is clear that the behavior of the critical frequencies depends on aspect ratio displaying characteristic functional forms for convective and rotational effects. From the analysis of Pr ¼ 0 it is seen that rotational effect becomes dominant as a is increased beyond 0:85. Further analysis is required for a41 to see if this trend continues over the aspect ratios. Three-dimensional time-dependent flow can be qualitatively reconstructed by summing the basic flow and the 3D

This research was supported by the Samuel and Anne Tolkowsky chair at the Technion-Israel Institute of Technology. The authors would like to acknowledge the use of computer resources belonging to the High Performance Computing Unit, a division of the Inter University Computing Center, which is a consortium formed by research universities in Israel. References [1] D.T.J. Hurle, Crystal Pulling from the Melt, Springer, Berlin, 1999. [2] S.M. Pimputkar, S. Ostrach, Convective effects in crystals grown from melt, J. Crystal Growth 55 (1981) 614. [3] Q. Xiao, J. Derby, Three-dimensional melt flows in Czochralski oxide growth: high-resolution, massively parallel, finite element computations, J. Crystal Growth 152 (1995) 169. [4] A. Yu Gelfgat, P.Z. Bar-Yoseph, A. Solan, Axisymmetry breaking instabilities of natural convection in a vertical bridgman growth configuration, J. Crystal Growth 220 (2000) 316. [5] J.C. Buell, I. Catton, The effect of wall conduction on the stability of a fluid in a right circular cylinder heated from below, ASME J. Heat Transfer 105 (1983) 255. [6] A. Yu Gelfgat, A. Rubinov, P.Z. Bar-Yoseph, A. Solan, Numerical study of three-dimensional instabilities in a hydrodynamic model of Czochralski growth, J. Crystal Growth 275 (2005). [7] W.E. Langlois, Buoyancy driven flows in crystal-growth melts, Ann. Rev. Fluid Mech. 17 (1985) 191. [8] A.A. Wheeler, Four test problems for the numerical simulation of flow in Czochralski crystal growth, J. Crystal Growth 102 (1991) 691. [9] A.T. Patera, A spectral element method for fluid dynamic: laminar flow in a channel expansion, J. Comput. Phys. 54 (1984) 468. [10] C. Cuvelier, C. Segal, A.A. van Steenhoven, Finite Elements Methods and Navier–Stokes Equations, D. Reidel Publishing Company, Dordrecht, 1986. [11] Y. Saad, M. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856. [12] G.H. Golub, C.F. van Loan, Matrix Computations, Johns Hopkins, Reading, MA, 1996. [13] R. Barrett, M. Berry, T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, SIAM Publications, Philadelphia, 1994. [14] A.Yu. Gelfgat, P.Z. Bar-Yoseph, A. Solan, T. Kowalewski, An axisymmetry-breaking instability in axially symmetric natural convection, I.J. Trans. Phenom. 1 (1999) 173.

ARTICLE IN PRESS Y. Rosenstein, P.Z. Bar-Yoseph / Journal of Crystal Growth 305 (2007) 185–191 [15] J.M. Lopez, F. Marques, J. Shen, An efficient spectral-projection method for the Navier–Stokes equations in cylindrical geometries II. Three dimensional cases, J. Comput. Phys. 176 (2002) 384. [16] M. Azaiez, A. Fikri, G. Labrosse, A unique grid spectral solver of the ND Cartesian unsteady stokes system, Illustrative numerical results, Finite Elem. Anal. Des. 16 (1994) 247. [17] R. Touihri, H. Ben Hadid, D. Henry, On the onset of convective instabilities in cylindrical cavities heated from below, i. Pure thermal case, Phys. Fluids 11 (1999) 2078. [18] J. Kevorkian, J.C. Cole, Perturbation Methods in Applied Mathematics, Springer, Berlin, 1981, p. 452.

191

[19] H. Tomonari, M. Iwamoto, K. Kakimoto, H. Ozoe, K. Suzuki, T. Fukuda, Standing-oscillatory natural convection computed for molten silicon in Czochralski configuration, Chem. Eng. J. 71 (1998) 192. [20] A. Yu Gelfgat, P.Z. Bar-Yoseph, A. Solan, Effect of axial magnetic field on three-dimensional instability of natural convection in a vertical bridgman growth configuration, J. Crystal Growth 230 (2001) 63. [21] A. Yu Gelfgat, P.Z. Bar-Yoseph, A. Solan, Three-dimensional instability of axisymmetric flow in a rotating lid-cylinder enclosure, J. Fluid Mech. 438 (2001) 363.