Use of a traveling magnetic field in VGF growth: Flow reversal and resulting dopant distribution

ARTICLE IN PRESS Journal of Crystal Growth 310 (2008) 3699– 3705

Contents lists available at ScienceDirect

Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro

Use of a traveling magnetic field in VGF growth: Flow reversal and resulting dopant distribution Ilma¯rs Grants , Gunter Gerbeth Forschungszentrum Dresden-Rossendorf, PO Box 510119, 01314 Dresden, Germany

a r t i c l e in fo

abstract

Article history: Received 22 November 2007 Received in revised form 14 May 2008 Accepted 16 May 2008 Communicated by A.G. Ostrogorsky Available online 24 May 2008

The melt flow in a model of a vertical gradient freeze facility under the action of a traveling magnetic field (TMF) is studied numerically. The radial temperature gradient over a concave solidification interface drives a radially converging flow which may cause an undesirable dopant concentration peak on the axis. We study the characteristics of such a flow for parameters beyond the previously reported linear regime. An upward directed TMF induces a body force which counteracts buoyancy. We report conditions under which the TMF reverses the flow direction for a wide parameter range. These conditions depend primarily on the product of the thermal gradient and the interface deflection. The simulation of the dopant transport demonstrates that the concentration peak disappears as soon as the flow direction changes. & 2008 Elsevier B.V. All rights reserved.

PACS: 44.27.+g 47.55.p 81.10.Fq 47.65.d Keywords: A1. Fluid flows A1. Segregation A2. Bridgman technique

1. Introduction Macrosegregation belongs to the most important problems concerning semiconductor crystal growth from the melt. In their recent review Hurle and Rudolph [1] described the Burton– Prim–Slichter (BPS) solution [2] as the single most useful piece of theory ever produced for practical crystal growers. No doubt, the BPS solution represented a milestone in the understanding of the concentration distribution in the melt and at the solidifying phase boundary of the growing crystal. The BPS approach, however, is only valid when the melt velocity is directed towards the solidification front [3] as it is typically the case for the crystal rotation driven flow in a Czochralski configuration. If the flow is in opposite direction, a lateral segregation unavoidably occurs causing a concentration peak along the axis of a converging melt flow [3,4]. Such a flow pattern takes place, e.g., over a concave solidification interface in vertical gradient freeze (VGF) growth [5] or if the flow is driven by external magnetic fields such as a rotating magnetic field (RMF) or a downward traveling magnetic field (TMF). In VGF growth, the  Corresponding author. Tel.: +49 351 260 3391; fax: +49 351 260 2007.

E-mail address: [email protected] (I. Grants). 0022-0248/$ - see front matter & 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2008.05.025

flow direction is reversed over a convex front and Hofmann et al. [6] indeed observed experimentally a much lower radial segregation in the latter case. The change from concave to convex front shape is only possible at the cost of a lower growth speed, provided that the furnace temperature distribution is fixed [7]. The flow direction can be reversed without changing the interface shape (and not loosing, thus, the productivity) by an upward TMF [8]. The main goal of our paper is to uncover the conditions of such a flow reversal caused by a TMF. To better understand the results we also perform a comprehensive numerical study of the natural buoyant flow. For a sufficiently small crystal radius the velocity of the natural flow scales linearly with the product of the thermal gradient and the solidification interface depth [5,9]. We extend these studies by considering the non-linear flow regime at higher parameter values. Besides, we also simulate the concentration field in order to illustrate the effect of flow reversal on the lateral segregation.

2. Problem formulation 2.1. Melt flow Considering the large number of independent parameters which may take place in different VGF configurations and trying

ARTICLE IN PRESS 3700

I. Grants, G. Gerbeth / Journal of Crystal Growth 310 (2008) 3699–3705

to keep the effort manageable and conclusions possibly general we followed the approach of Motakef [5,9]. Namely, we assumed a given shape of the solidification interface, a constant temperature gradient along the side wall and axial symmetry. Concerning the TMF we applied the most simple force expression [10] valid under low frequency and long wave conditions. The stationary Navier–Stokes equation in Boussinesq approximation for the velocity field Vðr; zÞ and pressure field pðr; zÞ writes ! =p soB2 a 2 2 ðVrÞV ¼ n= V  þ bgðT  T sol Þ þ r ez (1) r 8r

with boundary conditions  qH HjSþG ¼ ¼ 0, qn SþG TjS ¼ 0;

TjG ¼ z  d.

(2)

and no-slip conditions

V z ¼ H þ 0:5r

(3)

The temperature field Tðr; zÞ is described by the energy equation ðVrÞT ¼ kr2 T

(4)

qH ; qr

V r ¼ 0:5r

qH qz

(11)

and the operators L and N are defined as L¼

VjSþG ¼ 0.

(10)

This problem is governed by five dimensionless parameters which are summarized in Table 1. The flow velocities are expressed from the reduced stream function Hðr; zÞ by the following expressions [10]:

with the incompressibility constraint =V ¼0

(9)

q2 q2 3q ; þ þ qz2 qr 2 r qr

N ¼ Vz

q q þ Vr . qz qr

(12)

The function Wðr; zÞ represents the azimuthal vorticity as W ¼ r 1 ð=  VÞf . We are looking for solutions of the above problem with two aims:

with the boundary conditions (5)

 to determine the main properties of the natural flow (F ¼ 0) as

The meaning of symbols in the above equations is as follows: n; r; b; s, k and T sol denote kinematic viscosity, density, volumetric thermal expansion coefficient, electric conductivity, temperature diffusivity and solidification temperature of the melt, respectively. The TMF properties are the amplitude of the axial magnetic flux density B, the angular frequency o and the axial wave number a. The imposed thermal gradient is denoted by T 0 and the gravity acceleration by g. The solidification interface is denoted by S and the rest of the melt domain boundary by G, respectively (see Fig. 1). The molten zone has radius Ro and maximum length Lo . Eqs. (1)–(5) are further transformed into a dimensionless stream function and vorticity formulation [10] by applying the distance, time and temperature scales Ro , R2o =n and T 0 Ro , respectively:

 to determine the balance forcing F ¼ F b ðGr; d; Z; PrÞ such that

TjG ¼ T sol þ T 0 ðz  dÞ.

TjS ¼ T sol ;

NðWÞ ¼ LðWÞ  NðTÞ ¼ Pr

1

Gr qT  F, r qr

(6)

2

r T,

(7)

LðHÞ þ 2W ¼ 0,

(8)

2Ro z Lo Γ

a function of Gr, d and Pr for fixed values of Z. the second order derivative of the normal velocity is eliminated by the TMF near the bottom of the solidification 2 interface dip, i.e. qqzV2z ð0; 0Þ ¼ 0. This condition describes the exact instance of flow reversal at the solidification front.

2.2. Dopant advection To illustrate the effect of the flow direction change on the dopant concentration field we considered the advection problem under dilute mixture conditions. Under these conditions the melt flow does not depend on the dopant concentration. Given the velocity field Hðr; zÞ as an input, the dopant concentration field Cðr; zÞ is described by the mass conservation law qC þ NðCÞ ¼ Sc1 r2 C, qt

(13)

where Sc ¼ n=D is the Schmidt number and D is the dopant diffusion coefficient. The mass conservation at the solidification interface implies  qC  ð1  kÞSc Rep C þ  ¼ 0, (14) qz SðtÞ where Rep ¼ up Ro =n is the Reynolds number associated with the pulling rate up , and k is the equilibrium partition coefficient. The melt boundary G is assumed impermeable for the dopant  qC  ¼ 0. (15) qn G

Table 1 The governing dimensionless parameters

δ 0 Σ :z = δ (r/Ro

)2 Fig. 1. Scheme of the model.

Parameter

ΔT T−Tsol

Symbol

Definition

Range

Grashof number

Gr

bgT 0 R4o =n2

½103 : 108 

Magnetic forcing parameter

F

Prandtl number Aspect ratio Front deflection ratio

Pr Z d

soB2o R5o a=4rn2 n=k Lo =Ro d=Ro

½10 : 107  [0.002:0.5] [0.3:4] [0.01:0.3]

ARTICLE IN PRESS I. Grants, G. Gerbeth / Journal of Crystal Growth 310 (2008) 3699–3705

3. Numerical method 3.1. Melt convection The flow problem (6)–(10) was discretized by a spectral collocation approach and solved by the Newton method. For this purpose the curvilinear melt domain ðr; zÞ was first transformed into a rectangular ½0 : 1  ½1 : 1 domain for ðy; xÞ by the coordinate transformation y ¼ r;

x¼1þ

2ðz  ZÞ Z  dr

2

.

(16)

3701

TjGþS ¼ z. Note that the solution is not a trivial function of the mapped coordinates ðy; xÞ and it was obtained within the machine accuracy. Additionally we tested the implementation of the convective and buoyancy terms. As the second test we made a comparison with the results by the program for a regular cylinder which have been extensively validated during our former studies [12]. Instead of Tjz¼0 ¼ 0 (that would drive no flow) we applied the boundary conditions Tjz¼0 ¼ dð1  r 2 Þ, which for increasingly small d provide an increasingly similar driving buoyant force distribution in the vicinity of S. For d ¼ 0:02 the difference between the velocity fields by the different programs was below 1%.

The derivatives in transformed coordinates are q qx q ¼ ; qz qz qx

q qx q q ¼ þ , qr qr qx qy

(17)

4.1. Natural flow

where qx=qz and qx=qr are expressed from (16) as qx 2 ¼ ; qz Z  dr 2

qx rdðx  1Þ ¼2 . 2 qr Z  dr

4. Results

(18)

The solution was represented by its values on the Gauss– Lobatto mesh [11]. The fast Fourier transform (FFT) was then applied to express it as a series of Chebyshev polynomials. In this form the solution can be easily differentiated [11] with highest order accuracy. The nodal values of the solution derivative were then obtained by a backward FFT. Applying repeatedly the above approach to express all derivative terms in Eqs. (6)–(8), a set of discrete equations was finally obtained. This set required a zero residual at the internal nodes of the Gauss–Lobatto mesh [11]. The non-linearity of Eqs. (6)–(8) was resolved by the Newton iteration. The solution converged to a value of the residual below 106 within (typically) 4–6 iterations. The resolution was varied between 13  25 for Grp104 and 25  49 for GrX107. 3.2. Dopant advection An additional simplification was applied for the calculation of the dopant advection. We transformed the computational domain by letting y0 ¼ r, x0 ¼ 2z=Z  1 instead of Eq. (16) and replacing Hðy0 ; x0 Þ by Hðy; xÞ in the advection equation (13). In other words, we solved Eq. (13) in a physically rectangular domain with the velocity field mapped from the curvilinear one. The front curvature is a principal source of the flow and, therefore, cannot be avoided in the flow simulation whatever small it is. Concerning the dopant advection far from the diffusion dominated regime, in turn, the role of the curvature is limited to a geometric correction, which vanishes for d51. Eq. (13) was integrated in time by the explicit second-order Adams–Bashforth scheme [11]. The velocity field Hðy; xÞ was periodically updated with the steady solution of Eqs. (6)–(8) for the corresponding aspect ratios Z n ¼ Z init  nRep dt, with n ¼ 0; 1; . . . during the simulated growth. The step of interface position update was Rep dt ¼ 0:005. The initial condition of the concentration field was Cjt¼0 ¼ 1. The resolution of the concentration field was 129  97 modes—much higher than for the velocity field since thin solutal boundary layers may form near the solidification interface as well as internal diffusion layers between the vortices. 3.3. Validation The validation of the flow calculation was performed in two steps. First, we verified that the mapping is properly implemented by solving a test problem with a trivial solution Tðr; zÞ ¼ z; V ¼ 0. For this purpose we changed the boundary conditions (10) to

Motakef [5] found that the velocity field scales linearly with Gr  d provided that this quantity is small enough. After reaching a certain value of about 104 at Pr ¼ 102 a deviation from this dependency was observed. To explore the non-linear regime that follows, we calculated the maximum axial velocity at the axis for various Gr 2 ½104 : 108  and d 2 ½0:02 : 0:3 in a Z ¼ 2 cylinder for Pr ¼ 0:0072 (germanium). We observed that initially overlying curves of different d separated in the non-linear regime and each followed a Gr 1=3 scaling (Fig. 2a). Simultaneously, the vertical position z of the maximum velocity decreased proportionally to Gr 1=7 (Fig. 2d). To find out the flow dependency on d in the nonlinear regime we plotted Gr 1=3 V max ðdÞ (Fig. 2b) and obtained an approximately linear scaling for dp0:1. The position of the 1=7 velocity maximum z scaled as d (Fig. 2e). Finally, Fig. 2c shows 1 the Prandtl number dependency of the quantity Gr 1=3 d V max ðPrÞ 2=3 with a Pr scaling in the non-linear regime. The maximum position zðPrÞ followed a Pr1=5 scaling. Thus, the data in Fig. 2 suggest the existence of a similarity solution replacing the linear regime by the form V z ðzÞ ¼ V o f ðz=zo Þ where V o ¼ dGr

1=3

Pr 2=3

and

zo ¼ ðd=GrÞ1=7 Pr 1=5

in a wide range of parameters for a variety of VGF processes. To demonstrate this we plotted the reduced velocity V z =V o as a function of the reduced coordinate z=zo for three cuts in the parameter space through the point Gr ¼ 106 , d ¼ 0:1, Pr ¼ 0:0072 in Fig. 3(a)–(c), respectively. 4.2. Force balance estimates To estimate the TMF strength needed to counterbalance the buoyancy, let us equalize the characteristic buoyant force rbgT 0 d with the characteristic value of the TMF force [10] f / soB2 aR2o

(19) msoR2o

¼ 1 which when applied to at low frequency conditions Eq. (19) yield f / B2 a=m. We arrive at qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (20) B ¼ c bgT 0 dmr=a, where the unknown coefficient c is a function of four parameters c ¼ cðGr; Pr; Z; dÞ. Condition (20) can be expressed in terms of dimensionless criteria as F b ¼ 0:25c2 dGr from which c defines as rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F b ðGr; Pr; Z; dÞ . (21) cðGr; Pr; Z; dÞ ¼ 2 d Gr We defined the force balance as a state when the second derivative of the axial velocity vanishes at z ¼ 0; r ¼ 0. Note that the first derivative of V z is zero at this point due to no-slip and

ARTICLE IN PRESS 3702

I. Grants, G. Gerbeth / Journal of Crystal Growth 310 (2008) 3699–3705

10

1000 x

10 x −2

d = 0.02 d = 0.05 d = 0.10 d = 0.20 d = 0.30

1/3

x

1

103

105 Gr d

1

8 d = 0.01 7 log (Gr) = 6 5 4

108

Gr = Gr = 107 Gr = 106 0.1 0.01

0.1 10

x

1

/3

Vmax/Gr1/3d

10

Vmax/Gr1/3

Vmax

100

107

0.1 d

0.1 0.001

1

Gr = 106 d = 0.1 0.3

0.01

0.1

10

10

Gr = 106 d = 0.1 0.3

x −1/

7

d = 0.02 d = 0.05 d = 0.10 d = 0.20 d = 0.30

x

Gr = 108 Gr = 107 Gr = 106 1 0.01

0.1 103

10

105 Gr d

107

x −1/

5

ζ (Gr/d)1/7

1/7

ζGr1/7

ζ

1

1

Pr

0.1 d

8 d = 0.01 7 log (Gr) = 6 5 4 1 0.001

1

0.01

0.1

1

Pr

Fig. 2. Dependency of the maximum axial velocity (a–c) and characteristic boundary layer depth (d–f) on the governing dimensionless parameters; Pr ¼ 0:0072 for a, b, d and e.

Gr = 105 Gr = 106 Gr = 107 Gr = 108

0.25 reduced velocity

0.3

0.3

0.3

0.2 0.15

d = 0.30 d = 0.20 d = 0.10 d = 0.05 d = 0.02

0.25 0.2 0.15

0.2 0.15

0.1

0.1

0.1

0.05

0.05

0.05

0

0

0

−0.05

−0.05

−0.05 0

1

2

3

4

5

6

Pr = 0.002 Pr = 0.005 Pr = 0.02 Pr = 0.1 Pr = 0.5

0.25

0

1 2 3 4 5 reduced axial coordinate

6

0

1

2

3

4

5

6

Fig. 3. The dependency of the reduced axial velocity V z Pr 2=3 =Gr1=3 d on the reduced coordinate zðGr=dÞ1=7 Pr 1=5 at the axis ðr ¼ 0Þ for: varying Gr at Pr ¼ 0:0072, d ¼ 0:1 (a); varying d at Gr ¼ 106 , Pr ¼ 0:0072 (b); varying Pr at Gr ¼ 106 , d ¼ 0:1 (c).

incompressibility. This definition means that the flow direction close to this point reverses if F4F b . The velocity of front advancement Rep is typically much less than the velocity in the bulk and it is neglected in the definition of F b .

4.3. Numerical results of balance conditions The numerical results are presented in Fig. 4 in terms of the parameter c. This parameter varied little over a wide range of Gr

including both the linear and non-linear regimes. A slightly higher balance forcing (a larger value of c) was observed in the non-linear regime with large Gr. The cðZÞ dependency was constant for large enough Z and rapidly decreased for small Z. The maximum of cðZÞ was observed at Z ¼ 1:65 or 1.1 for Gr ¼ 106 or 107, respectively (Pr ¼ 0:0072, d ¼ 0:1). For a higher Prandtl number of 0.026 and d ¼ 0:1 the position of the maximum was Z ¼ 1:07 or 0:753 at Gr ¼ 106 or 107, respectively. The case of Gr ¼ 103 shown with triangles in Fig. 4b represents the linear flow regime. In this regime c did not depend noticeably on Gr and Pr.

ARTICLE IN PRESS I. Grants, G. Gerbeth / Journal of Crystal Growth 310 (2008) 3699–3705

2.8 d = 0.30 d = 0.10 d = 0.02

2.6

3703

0.75

0.75

0.5

0.5

0.25

0.25

2.4 0

0.5

1

2

0

0.5

1

2

2

1.5

1.5

1

1

0.5

0.5

height

c

0 0

2.2

1.8 1.6 2.5

3

3.5

4

4.5

5 5.5 log (Gr)

6

6.5

7

7.5

3 2.8 2.6 0

c

0

0.5

1

2.4

0 0

0.5

1

radius

Gr = 103 Gr = 106 Gr = 107 Pr = 0.0072 Pr = 0.026

2.2 2

Fig. 5. Isolines of the stream function [0:5rHðr; zÞ] for the natural flow (a and c) and the corresponding TMF balanced state (b and d) for Z ¼ 0:75 (a and b) and Z ¼ 2 (c and d); Gr ¼ 107 , Pr ¼ 0:026; The streamline step is 2.

1.8 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

TMF strength increased (Fig. 6b, Fp105 ). The peak was eliminated as soon as the flow direction changed.

Z Fig. 4. The balance conditions for: varying Gr and fixed Z ¼ 2, Pr ¼ 0:0072 (a); varying Z and fixed d ¼ 0:1 (b).

The effect of a TMF with F ¼ F b on the resulting flow field is illustrated in Fig. 5. Two cases are shown with Z ¼ 0:75 and 2 which correspond to the maximum of the cðZÞ dependency and a ‘‘large’’ Z value, respectively.

4.4. Dopant concentration Fig. 6 illustrates the effect of a TMF of various strength on the dopant concentration in the crystal. The parameters corresponded to a gallium doped germanium growth (k ¼ 0:09, Sc ¼ 20) at a dimensionless pulling rate of Rep ¼ 0:5. The other parameter are fixed to Gr ¼ 106 , d ¼ 0:1 and an initial height of Z init ¼ 3. One of the cases considered (F ¼ 1:5  105 ) was very close to F b ðZ init Þ ¼ 1:49  105 and initially there is F4F b , that means a flow towards the front. As for ongoing crystal growth the aspect ratio of the liquid zone approached the maximum of the cðZÞ curve (Fig. 4b), the flow direction turned away from the solidification interface. This event is marked in the axial concentration profile by a pronounced spike. The radial concentration profiles showed a pronounced concentration peak for all cases of a converging (i.e., buoyancy dominated) flow. The height of this peak initially increased as the

5. Discussion We observed two regimes of the natural low-Pr flow above a concave solidification interface: The linear regime with a Pr independent velocity of V max  4:1  103 Gr  d and the nonlinear regime with V max  0:25Gr 1=3 dPr 2=3

at z ¼ z  1:8ðd=GrÞ1=7 Pr 1=5 .

(22)

A self-similar asymptotic solution seems to exist in the latter regime. Let us estimate the border between those regimes by equalizing both numerically observed scalings of V max . That yields Gr max  400=Pr

or

Ramax  400,

(23)

independently on d, where the Rayleigh number is Ra ¼ GrPr. For instance, the linear regime persists up to Gr max  5  104 for gallium that translates into Ro  0:7 cm under a 3 K/cm thermal gradient. Thus, the linear regime may be expected for small crystals only. Condition (23) can be also expressed as Pe  1:6d, where Pe ¼ Pr V max is the Peclet number. The latter condition puts the linearity border (23) at an increasingly low velocity as d is decreased (the effect is seen in Fig. 2a). Thus, this regime change cannot be due to the non-linear convection term NðWÞ in Eq. (6) as commonly believed. The form of (23) suggests the braking force of buoyancy, instead. The melt convection deforms the temperature field by an amount measured by the Peclet number. This thermal field deformation creates a buoyancy force which tries to eliminate the flow. The regime change takes place as soon as the

ARTICLE IN PRESS 3704

I. Grants, G. Gerbeth / Journal of Crystal Growth 310 (2008) 3699–3705

0.8

1.4 F=0 F = 1.5x105 F = 3.0x105 F = 0.5x105

1.2 concentration

1

0.8 F=0 F = 0.5x105 F= 1.0x105 F = 2.0x105 F = 3.0x105

0.7 0.6

0.8

z/Zinit = 0.10 = 0.50 = 0.75

0.7 0.6 0.5

0.5 0.6

0.4 0.4

0.4

0.3

0.3

0.2 0

0.2

0.2 0

0.2

0.4 0.6 z/Zinit

0.8

1

0.1 0

0.2

0.4 0.6 radius

0.8

1

0

0.2

0.4 0.6 radius

0.8

1

1

1

0.5

0.5

0 0

0.5

1

heigth

Fig. 6. The axial (a) and radial (b and c) concentration profiles for Gr ¼ 106, Pr ¼ 0:0072, d ¼ 0:10: r ¼ 0 (a); z=Z init ¼ 0:5 (b); F ¼ 1:5  106 (c).

0 0

0.5

1

radius Fig. 7. Isolines of the buoyant vorticity source ð1=GrdÞð1=rÞðqT=qrÞ in the non-linear (a) and in the linear (b) regime at Gr ¼ 107 and 103 , respectively. Isoline step 0.2, negative isolines dashed, the zero isoline bold.

braking convection effect on the thermal field (measured by Pe) becomes comparable to the flow-driving imposed radial gradient (measured by d). The buoyancy, thus, has a dual—driving and braking—role in the non-linear regime. This is well demonstrated in Fig. 7 where the buoyant vorticity source is compared in both regimes. The numerically observed velocity scaling held over a large range of the Grashof numbers. Some deviation from the scalings occurred when the position of the velocity maximum z became comparable to the interface depth d (circles in Fig. 2a and d at large Gr). In that case the scaling overestimated the velocity and underestimated z (see also Fig. 2b and e). The strongest impact on a deviation from scalings (22) had the dimensionless front interface depth d in the range of relevance for crystal growth. A larger deviation was observed for a larger Gr (Fig. 2b and e). These observations can be explained by the geometric correction becoming increasingly important when the characteristic thickness of the convection torus gets comparable to the dip depth. As the Prandtl number decreases the flow velocity increases and at a some point the convective terms should become important in Eq. (6). To estimate the border of the non-linear similarity solution regime, let us estimate the convective term NðWÞ  V 2max =z3 and equalize it to the driving buoyant force dGr. By making use of the numerical results (22) after some algebra we obtained the estimate of the minimum Prandtl number 0:78

Prmin  0:002Gr 0:13 d

.

(24)

According to this estimate Prmin ¼ 0:002 for Gr ¼ 106 and d ¼ 0:1. A considerable deviation from the similarity solution is indeed seen in Fig. 3c for this parameter combination while a higher Pr value provides a much better agreement. A characteristic Prandtl number in VGF is about 0.01 while the interface deflection dp0:25. Thus, the similarity solution should be valid in VGF for up to Gr ¼ 109 at least, according to the above estimates. The main focus of our study was on the conditions of equilibrium of natural and magnetic driving forces near the phase interface. It turned out that these conditions are surprisingly well approximated by the simple expression (20) derived from order of magnitude considerations. Namely, the shape function cðGr; Pr; Z; dÞ varied little over a wide range of the governing parameters. The flow regime change was observable also in the cðGrÞ dependency (Fig. 4). According to Eq. (23) the regime change took place independently on d at Gr  5  104 for Pr ¼ 0:0072. This corresponds to the minimum of the lower curve in Fig. 4. As seen in this figure, a relatively steep rise began at that point for all three cases presented. Pronounced maxima were observed in the cðZÞ dependency (Fig. 4b). These maxima Z max correlated with the characteristic length-scale so that Z max =zo 2 ½5: 6 for the four displayed cases. As seen in Fig. 5a, Z max corresponds approximately to 1.5 thicknesses of the natural convection cell. Let us consider the example of a 2 in germanium crystal with a 3 K/cm thermal gradient in the melt and d ¼ 0:1. With s ¼ 1:6  106 S=m the low frequency condition msoR2o ¼ 1 constrains the TMF frequency at o  800 s1 (125 Hz). The Grashof number is about 6:43  106 for which c  2:6 (b ¼ 1:1  104 K1 , n ¼ 1:35  107 m2 =s). Inserting now the germanium properties in Eq. (20) and assuming for instance a ¼ 40 m1 we have B ¼ 0:94 mT or F b ¼ 1:1  106 (r ¼ 5500 kg=m3 ). Thus, the required strength of the TMF is rather moderate if the flow direction change is the purpose. The flow reversal efficiently removed the concentration peak at the axis. If, in turn, the TMF strength was insufficient to change the flow direction (FoF b  1:5  105 ), the peak was enhanced (Fig. 6b). In the latter case the TMF just slowed the natural flow near the solidification interface. A weaker flow means less mixing, thus causing a larger variation of concentration. The radial variation of concentration decreased sharply as soon as the flow direction was reversed despite the velocity magnitude near the front was small as well. Note that the reversed flow eliminated the lateral segregation only in the central part. The upward velocity occurred now at the rim and is accompanied by the lateral segregation there. In this respect one may ask: Why the variation of concentration was so different in the two cases of F ¼ 105 and

ARTICLE IN PRESS I. Grants, G. Gerbeth / Journal of Crystal Growth 310 (2008) 3699–3705

If the knowledge of the thermal conditions is very approximate, then a TMF of variable strength can be used as a measurement/adjustment tool. To demonstrate this we performed a numerical simulation of segregation under the TMF with a gradually decreasing strength. The concentration peak started to grow causing a steep rise in the axial profile as soon as the buoyancy took over the TMF. These positions are marked by dotted lines in Fig. 8. Knowing the strength of the TMF at the beginning of the concentration peak one can recalculate the buoyancy force which is proportional to the product of interface deflection and the thermal gradient.

1.4 Finit =

1.2

concentration

2x105

Finit = 4x105 Finit = 8x105

1

3705

0.8 0.6 0.4

6. Summary

0.2 0 0

0.5

1

1.5 height

2

2.5

3

Fig. 8. The axial concentration profiles for a linearly decreasing in time TMF with an initial value F ¼ F init at z ¼ 0 and vanishing F at z ¼ 3; Gr ¼ 106 , d ¼ 0:1, Pr ¼ 0:0072, Z init ¼ 4. The dotted lines mark positions where F ¼ F b .

2  105 ? The main difference lies in the extent of the front portion with an upward velocity where the BPS approach fails [3]. The vertical velocity changes its sign at about r  0:8 (Fig. 5b and d). Thus, with FoF b the rejected dopant accumulate from r  0:8 to 0, while for F4F b it happens over a much shorter distance from r  0:8 to 1. Therefore, the radial segregation is much lower for F4F b and it is limited to the rim of the crystal. The above effect of the TMF on the radial segregation with F4F b was similar to the effect of the flow direction change by a convex solidification front [6]. The advantage of the TMF is that more latent heat is allowed by the concave solidification front. Concerning the axial segregation the velocity should be low to keep the effective segregation coefficient possibly close to unity. It may be a very challenging task since it requires the precise knowledge and control of the thermal conditions. Besides, such a fine balance imposes strong requirements for the rotational symmetry of both the thermal conditions and the TMF force. The axisymmetric parts of buoyancy and TMF forces eliminate each other when F  F b , but any uncontrollable asymmetry stays the same. Thus, the role of asymmetry increases as the forcing approaches the balance conditions leading likely to an essentially uncontrollable three-dimensional flow pattern.

The axisymmetric flow in a simplified model of VGF is studied numerically. The buoyancy flow scalings (22) with respect to Grashof number, interface deflection and Prandtl number are found for intermediate parameter values bounded by expressions (23) and (24). This regime is expected when the crystal diameter is roughly between 15 and 150 mm. The non-linear dependency of the velocity vs. the temperature gradient is caused by the braking buoyant force of the convection deformed temperature field. A traveling magnetic field with strength (20) offsets the natural flow near the center of the solidification interface. The coefficient c in this expression depends little on the governing parameters and its value lies between 1.5 and 3. The flow reversal removes the concentration peak in the central part and reduces the radial segregation significantly. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]

D.T.J. Hurle, P. Rudolph, J. Crystal Growth 264 (2004) 550. J.A. Burton, R.C. Prim, W.P. Slichter, J. Chem. Phys. 21 (1953) 1987. J. Priede, G. Gerbeth, J. Crystal Growth 285 (2005) 261. J. Priede, G. Gerbeth, Int. J. Heat Mass Transfer 50 (2007) 216. S. Motakef, J. Crystal Growth 102 (1990) 197. D. Hofmann, Th. Jung, G. Mu¨ller, J. Crystal Growth 128 (1993) 213. I. Grants, G. Gerbeth, J. Crystal Growth 207 (1999) 138. S. Yesilyurt, S. Motakef, R. Grugel, K. Mazuruk, J. Crystal Growth 263 (2004) 80. S. Kaddeche, J.P. Garandet, C. Barat, H. Ben Hadid, D. Henry, J. Crystal Growth 158 (1996) 144. I. Grants, G. Gerbeth, J. Crystal Growth 269 (2004) 630. C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral Methods in Fluid Dynamics, Springer, Berlin, 1998. I. Grants, G. Gerbeth, Phys. Fluids 16 (2004) 2088.