Convection and segregation in directional solidification of dilute and non-dilute binary alloys: Effects of ampoule and furnace design

Journal of Crystal Growth 80 (1987) 155—190 North-Holland, Amsterdam

155

CONVECTION AND SEGREGATION IN DIRECTIONAL SOLIDIFICATION OF DILUTE AND NON-DILUTE BINARY ALLOYS: EFFECTS OF AMPOULE AND FURNACE DESIGN Peter M. ADORNATO and Robert A. BROWN Department of Chemical Engineering and Materials Processing Center, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA

Received 18 April 1986; manuscript received in final form 7 July 1986

The effect of furnace configuration and ampoule design on the temperature field, the convection in the melt, the shape of the melt-solid interface, and the segregation of solute in the crystal are analyzed for the directional solidification of several dilute and non-dilute binary semiconductor alloys. The analysis is based on numerical calculations using a Petrov—Galerkin/finite-element method for solving the free-boundary problem that describes axisymmetric convection in the melt, the interface shape, and heat transfer in melt, crystal, and ampoule in a quasi-steady state model of the vertically stabilized Bridgman-Stockbarger system and for a furnace with a uniform temperature gradient imposed over a long section of ampoule. The flow in molten germanium grown in the Bridgman-Stockbarger system has two vertically-stacked toroidal cells. The top cell moves melt upward along the ampoule wall and is driven by the radial temperature gradients caused by the junction of the adiabatic and hot zones of the furnace. The cell adjacent to the melt/solid interface convects melt upward along the axis of the ampoule and exists because of the mismatch between the thermal conductivities of melt, crystal and ampoule. The relative intensities of the two cells depend on the values of the three conductivities and the magnitude of the axial temperature gradient. The upper cell is not present in the furnace with a uniform axial gradient; however, radial temperature gradients near the melt/crystal interface still drive convection. Adding a solute that is less dense than the bulk melt and that is preferentially incorporated upon solidification, such as silicon in germanium, decreases the intensity of the flow near the interface by increasing the melt density there. Radial temperature gradients interact with the solute field to cause the sideways diffusive instability predicted by Hart [J. Fluid Mech. 49 (1971) 279). Calculation of the melt/crystal interface shape, radial segregation, and effective segregation coefficient for the growth of dilute gallium is germanium agree reasonably well with measurements of Wang [PhD Thesis (1984)] when an accurate description of heat transfer in the ampoule and furnace is included in the analysis.

1. Introduction Interest in the details of heat and mass transport during directional solidification of semiconductors inside an ampoule has intensified in the last several years in an effort to quantify crystal growth in this geometry. Most theoretical research has focussed on understanding heat transfer between melt, crystal, and ampoule and has led to several important discoveries necessary for steadystate crystal growth. One-dimensional thermal analyses [1—4]have established the ampoule length necessary so that transients introduced by thermal end effects can be neglected. Two-dimensional heat transfer calculations have explored the effects of the ampoule on the curvature of the melt/solid interface [5—8j.The influence of subtle changes in

the heat transfer on the convection in the melt and its effect on solute segregation in the crystal has received much less attention, but is extremely important because it is the compositional uniformity, along with the crystallographic quality, that sets the value of the crystal. The purpose of this paper is to present detailed calculations coupling heat transfer in the melt, crystal and ampoule to convection in the melt and segregation in the growing crystal. We consider a binary alloy and account for convection driven by both temperature and concentration gradients. The temperature variations important in directional solidification are the axial gradient imposed for solidification and the radial gradients caused by imperfections in the heat transfer between the furnace, ampoule, melt and crystal. Thermal con-

0022-0248/87/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

156

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification

vection is caused by the radial gradients and by axial gradients that lead to an unstably stratified melt, Convection driven by concentration-dependent density differences is initiated in the diffusion ]ayer adjacent to the melt/solid interface caused by solute partitioning between these phases. In non-dilute alloy systems the interaction of the solutal and thermal driving forces is so closely coupled through convection of the solute field, that only fo!.~ numerical simulation can give quantitative insight into the structure of the flow, Previous calculations of melt flow in directional solidification have been limited to models with idealized boundary conditions that poorly model the temperature gradients in an actual furnace [9—13].Chang and Brown [10] and McFadden et al. [19] computed flows in an ampoule resulting from perfectly vertical temperature gradients. The first analysis concentrated on flows driven by an unstable temperature gradient and the latter research considered flows driven by the solute field in the presence of a stabilizing axial temperature gradient. In both cases flows evolved from a static state above critical values of the appropriate Rayleigh numbers, measuring the driving forces for convection. Introducing into the model calculations the radial temperature gradients present in real crystal growth systems leads to convection for all values of Rayleigh number and changes the structure of flows driven by the unstable axial gradients [10]. Even in the presence of stabilizing axial gradients, i.e. melt above crystal in directional solidification, the flows induced by radial temperature gradients can be extremely intense. Chang and Brown [9,11] and Carison et al. [12] examined this situation for idealized Bridgman—Stockbarger systems where thermal boundary conditions were imposed directly on the melt and crystal without any accounting for the presence of the ampoule. Both sets of calculations demonstrated intense thermal convection driven by changes in the temperature distribution along the ampoule boundary. In practice, the ampoule and the imperfect heat transport between the furnace wall and the ampoule surface distort the temperature profile imposed by the furnace. The axial gradients are

decreased by the presence of the ampoule and the radial gradients can be either increased or decreased, depending on the details of the heat transport between the charge and ampoule. Accurate calculations of convection and the concomitant solute segregation can only be modelled when the ampoule is included. Incorporating the ampoule into models like those reported in refs. [11,12] also makes possible the comparison of the predictions of calculations directly with experimental measurements in well characterized growth systems. Optimization of the ampoule and furnace design to decrease radial and axial solute segregation is also feasible. These are the goals of our research. Experiments on the growth of semiconductor alloys in small-scale vertical Bridgman—Stockbarger furnaces have advanced considerably using the existing theoretical base to optimize heat transfer. Two furnace designs are of particular interest in this study and are shown schematically in fig. 1. These are the classical Bridgman—Stockbarger furnace constructed by Wang and Witt [14] at the Massachusetts Institute of Technology and the furnace constructed by Rouzaud and Favier [15] at the Centre d’Etudes Nucléaires de Grenoble. The MIT furnace is based on heat pipes for forming isothermal hot and cold regions separated by a nearly adiabatic zone. The isothermal zones are long enough that the crystal growth rate equilibrates with the ampoule displacement rate when a dilute alloy is solidified. This condition is not satisfied in conventional furnaces where the crystal is extracted directly into the ambient [16]. The Grenoble furnace uses a tapered heating element to establish a nearly linear temperature profile over a length of ampoule approximately 30 times the ampoule radius. Again, for the growth of a dilute alloy, the crystal growth rate reaches the ampoule displacement rate after an initial thermal transient. When the alloys is not dilute the melting point of the material may vary appreciably with changes in the average solute concentration inside the ampoule, and lead to transients in the crystal growth rate. This transient has been calculated by Bourret et al. [17] in the limit of diffusion-controlled solute transport in the melt. Convective mixing in the melt will cause a continuous change

P.M. A dornato, R.A. Brown

/

Convection and segregation in directional solidification

Ampoule

157

Ampoule TH

w

o NJ

I II— -J

0

I

I

w =

\\\

\\•\

~L9

0

_____

t

L

00

0

j

\\\ __________

-J U) C 0 NJ

U

Cl)

a

_____

T~b

Fig. 1. Schematic diagrams of the models for the MIT (a) and Grenoble (b) vertical Bridgman crystal growth systems.

in concentration at the interface and a corresponding change in the melting temperature of a non-dilute alloy, Comparison between the predictions of calculations, like those reported here, and small-scale growth experiments is only possible when the thermal boundary conditions for the experiment are established precisely and when solute segregation is predicted by the calculations. The experiments described in refs. [14] and [15] meet the first criterion. The calculations of solute transport in these crystal growth systems are extremely difficult because the small solute diffusivity in the melt results in convection-dominated transport, even at low flow velocities, so that the Peclet number for species transport is large. Solution of

the species transport equation for high Peclet numbers is one of the outstanding problems in numerical analysis [18,19]. In non-dilute alloys, the solute is coupled to the flow through the concentration dependence of the buoyancy force. Then difficulties with numerical solution of convection-dominated transport also limit the calculation of the velocity field. We have extended the finite-element/Newtoniteration method development in refs. [9] and [10] to include a Petrov—Galerkin formulation of the field equations which is designed to aid in solution of convection-dominated transport equations. The description of the Petrov—Galerkin formulation and tests of the accuracy of the method for solving flow problems in directional solidification are in-

158

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification

cluded in ref. [20].The method is described briefly in section 3. Results are presented in sections 5 and 6 for the growth of both gallium-doped germanium and silicon—germanium alloys in the furnace configurations developed at MIT and Grenoble. The growth of gallium-doped germanium has been studied experimentally in both systems [14,15], and the group at Grenoble also has analyzed the growth of SiGe alloys. Calculations presented in section 5 for growth of gallium-doped germanium in ampoules of boron nitride, quartz, graphite, and a quartz/graphite composite ampoule show the effects of ampoule conductivity on flow structure and intensity. The impact of these changes on the solute concentration in the crystal is demonstrated by radial and axial segregation calculations for the dilute gallium—germanium alloy. The predicted values of radial segregation and effective segregation coefficient for gallium are compared to values obtained by Wang and Witt [14] in MIT. These are the first quantitative comparisons with experiments of solute segregation predicted from detailed analysis of melt flow, Adding silicon to germanium melt causes a stabilizing axial density gradient adjacent to the melt/crystal interface, as recognized by Rouzaud et al. [15] and damps convection driven by radial temperature gradients. This effect is demonstrated in section 6 by calculations for the constant gradient furnace of the Grenoble group.

2. Pseudo steady-state model The calculations described here are based on a pseudo steady-state model (PSSM) of a directional solidification system in which the translation of the ampoule is modelled by supplying melt at the top of the ampoule with a uniform velocity Vm and composition c0 and removing crystal from the bottom of the ampoule at the same average cornposition and growth rate Vm(pm/Ps), where (Pm, Vm) and (~V~)are the densities and velocities of the melt and crystal, respectively, Employing the PSSM neglects transients in the field variables (velocity, pressure, temperature, and concentration) caused by the steady decrease of =

=

the length of the melt with crystal growth and the displacement of the ampoule in the furnace. For sufficiently long ampoules the thermal end effects are negligibly small [2—4]. Axial segregation of solutes in directional solidification is due to the continual evolution of the amount and distribution of solute in the melt with length solidified. The changes in concentration are usually measured from the axial concentration profile along the crystal and reported as the evolution of the effective segregation coefficient keff with distance grown [16]. The composition of the solid changes continuously when strong mixing is assumed in the melt, as accounted for in the well-mixed model of Scheil [211. In the limit of a well-mixed melt, the effective segregation coeff icient is defined in terms of the concentration field in the melt as k

—

k


—

—

argue that the first stage of solute redistribution is rapid convection followed by the much slower diffusion process. For systems with the intense convection and steep concentration gradients characteristic of our previous calculations for Bridgman growth [9], convective mixing occurs quickly and diffusion is needed only over very small boundary layers to smooth the concentra-

P.M. Adornato, R.A. Brown

/ Convection and segregation

159

in directional solidification

tion profile. Then the PSSM is expected to be

where hr

and (er, e0, e~) are the unit vectors in cylindrical coordinates, Convection in a binary melt driven by density gradients due to both temperature and composition is most simply described by the linear relationship

valid. The experiments of Wang and Witt [14] confirmed this assumption, as discussed below. The PSSM is equally valid in the limit of diffusion-controlled growth in melts long enough that the diffusion layer adjacent to the melt/crystal interface does not reach the end of the ampoule. One of the purposes of crystal growth in microgravity has been to achieve these conditions. The relevance our segregation calculationsintoreal thecrystal prediction of axial and of radial growth systems is discussed in section 7 in terms of the convection-controlled and diffusion-controlled limits.

~ T~0t)+ $~(c c0)], (4) where the reference density is that of melt at the Thot and composition c 3T and temperature 0, and / f~are the coefficients of thermal and solutal expansion, respectively. Eq. (4) is written in dimensionless form as as p 1 Pr [Rag S RaT 9] (a~/gL3),

dh/dr

—

—

—

=

—

(5) 2,1. Equations and boundary conditions

where p is the dimensionless density (~/~) and Pr, RaT, and Ra~are the Prandtl number, and

The field equations and interface shape are described in terms of the stationary cylindrical polar coordinate system (r, z, 0) with origin located at the centerline and the top of the ampoule. Variables are put in dimensionless form by scaling lengths with the length of the ampoule L, velocity v(r, z) with the thermal diffusion velocity of the melt am/L, pressure p(r, z) with Pm~n/1~2 and composition with the inlet value c

thermal and solutal Rayleigh numbers. These and other dimensionless groups are defined in table 2. The relative effects of temperature and cornposition on the bulk motion of the melt depend on the direction of the solute gradient and the signs of the Rayleigh numbers. In our convention RaT is positive if the melt is positioned above the crystal relative to a gravitational vector —ge’.

0.

Then the radially-averaged density decreases with increasing vertical distance and is stabilizing. When viewed as a perturbation from diffusioncontrolled growth the primary solute gradient that drives flow is due to the axial partitioning of

Dimensionless concentration and temperature fields are defined as S(r, z) as c(r, z)

e(r

z)

= —

—

1,

T( r, z) — T cold (2) ThOt—l~O1d 1~oId are reference temperatures

where and design under consideration. The for the ThOr furnace set point temperature of the hot and cold zone are used for these values in the MIT furnace; the temperatures at the ends of the constant gradient are used in the Grenoble furnace, The location of the melt/crystal interface is described by z h(r), 0 r ~ A, where A as R/L is the aspect ratio of the ampoule. The unit vectors normal N and tangential t to the surface are written as =

N~ ___________ ez~hrer (i + h~)1~2 ~

~

er+hrez + h~)’~2

___________

(i

(3)

solute between melt and crystal. Whether or not this concentration gradient stabilizes the density profile depends on the density of the solute relative to that of the bulk melt and on the direction of the partitioning, i.e. whether k is less than or greater than unity. We consider the crystal growth of dilute gallium—germanium and of non-dilute silicon—germanium alloys. Solutal convection is neglected in the first system. Silicon is the lighter component (/3~< 0, Ra <0) and partitions into the crystal in SiGe, thus leaving heavier melt adjacent to the interface, so that the dependence of the density on the silicon concentration is expected to stabilize convection near the interface. The axisymmetric forms of the Boussinesq equations are (6)

160

P.M. Adornato, R.A. Brown

v~V v

=

—

Vp + Pr

28, v~vO v (Sc/Pr) v. VS =

=

v

—

Pr (Ra5 S

—

/ Convection and segregation

RaT 0) e~, (7) (8)

v 2~

(9)

where V e8/8r + e~8/azand the Schmidt number Sc is defined in table 2. The energy equations in the crystal and ampoule are Pec ~ v0 ~ 2~ (10) =

Pea e~’V9

(11)

‘YaV20,

=

where ‘/c ac/am and ‘Ya as aa/am are ratios of the thermal diffusivities of each phase to the value for the melt. The Peclet number Pe is the dimensionless transition rate of the ampoule scaled with the characteristic velocity for thermal diffusion; see table 2. The shape of the melt/crystal interface is set by the slope of the liquidus curve em(S) from the binary phase diagram as 9(r, h(r))

Urn(S) as 9~+ m(S + 1

=

—

1/k), (12)

where the curve is approximated by a straight line with dimensionless slope m and 0,~is the melting temperature of the alloy with concentration c 0. We assume that m and the segregation coefficient k are independent of concentration; these are good assumptions for SiGe at low weight percent Si. The numerical method presented below does not require these simplifications. The interfacial temperature and solute fields are related to the corresponding values in the bulk melt by the balances [N. V0]m — K~LN’vU],.

[N.vS]=

=

St Pe (N. e.),

Pe Sc (e~’N)(1—k)(S+1), Pr

‘

t =

Pe (e~ t), ‘

G(v N)

=

Pe (e, N), .

that v~ 0 and v~ Pe there. When ~ is not unity, the second condition in eq. (14) dictates an aphysical in the velocity atto the ampoule wall. This jump condition corresponds the presence =

=

of a singularity at the contact curve between melt, crystal, and ampoule. Analysis of a similar singularity for liquid/liquid/solid contact lines has demonstrated the necessity of allowing the fluid to slip in a small region near the contact line [23,24]. Because we deal with materials having nearly equal melt and solid densities, we avoid introducing melt slip by setting the densities to be equal, i.e. = 1. The thermal boundary conditions shown in fig. 1 are specified along the ampoule wall according to the heat transfer condition [aO/ar]a

=

Bi(z) [9— 9~(z)],

(16)

where Bi( z) is a dimensionless heat transfer coefficient Bi haL/ka defined to include radiative, conductive, and convective transport between the ampoule and the furnace, and 0~( z) is the furnace temperature distribution. Applying the form of eq. (16) shows above assumes that the ampoule and furnace temperatures are not very different at a specific axial location so that the radiative contribution can be linearized, as was done by Jasinski

(13)

et al. [4]. The particular profiles for Bi(z) and Oas(Z) used to model each furnace are discussed in section 5. The specification of the solute concentration field is completed by setting the diffusive flux of solute through the inlet to be equal to PeS so that the PSSM model is equivalent to the transport of bulk melt with the dimensional average concentration of c0. Perfect thermal contact at the melt/ampoule and crystal/ampoule interfaces is assumed, so that the conductive fluxes at these boundaries are

(14)

equal:

where K,. as kc/km is the ratio of thermal conductivities between crystal and melt. The boundary conditions on the velocity field at the interface are o

in directional solidification

(15)

where a p,./p~. The top and sidewall of the ampoule are assumed to be no-slip surfaces so

~ !.

3r

r

=

.i

A.

=

Ka[~~,

K

= c~

Ka{~j~

arj (17)

The axis of the cylinder is taken as a line of symmetry for all field variables, as is consistent with the assumption of axisymmetric fields. The top and bottom of the ampoule are assumed to be

P.M. Adornato, R.A. Brown

/

161

Convection and segregation in directional solidification

at the temperature of the hot and cold zones, i.e. 0 = 1 at z = 0 and 9 = 0 at z = 1. These conditions have no influence on the calculations when the ampoule is sufficiently long. The velocity fields computed from the solution set eqs. (6)—(17), are displayed as contours of the

temperature gradient to change directions over short length scales in the ampoule. The magnitude of the radial temperature gradient must vary proportionally with the vertical gradient defined for the Bridgman system, i.e. L~T*/R LXT/L5 where z~T=ThOr ~oId’ Using

stream function ‘I’(r, z) defined by the relations

this definition in eq. (18) relates Gr to the Rayleigh number as 3, RaT=GrPr/A

VZ

ia’i’ r8r’

We use these definitions in section 2.2 to define relative measures for the intensity of the convection. 2.2. Scaling analysis

When melt is positioned above the ampoule flows are driven in either geometry solely by the radial temperature gradients in the melt. The radial gradients drive cellular motions, which are superimposed on the uniform flow caused by crystal motion. When the thermal convection is strong enough to result in melt recirculation, the relative strength of this flow to the uniaxial convection can be computed by comparing the magnitude of the stream function at the center of the recirculation ~1~*(r*, z*) to the value of ~I’at the wall ‘Ic. The ratio (~I~*/,IF measures the volumetric flow rate of melt in the recirculation relative to the flow rate needed for crystal growth at the speed Vm across the vertical cross section )

Z = Z ~.

—

—

where Gr = g$T(~T/Lg)R/v.

(19b)

The Grashof and Rayleigh number appropriate for directional solidification in a stabilized axial gradient are proportional for fixed values of the length of the gradient zone, the ampoule length and the thermophysical properties. For an arbitrary furnace and ampoule design, the convection pattern in the melt will depend on both length scales L and R, so the choice of scaling for the convective driving force can be described equally well in terms of either Gr or RaT. For the well-designed furnaces of the MIT and Grenoble groups, the flow structure is independent of the melt length beyond a minimum value, so that Gr is the most appropriate measure. Since the ampoule length is held constant in our calculations RaT and Gr can be used interchangeably. Scaling analysis is useful for elucidating some

The the 1inai Rayleigh number RaT defined in terms of th... vertical temperature difference and the length of the ampoule does not scale directly the radial gradients and hence this dimensionless group cannot be used directly in order-to-magnitude estimates of the driving forces for recirculating convection. The true driving force should be measured on the basis of a characteristic radial temperature differences ~ T * across the melt in terms of a Grashof number 4/v2. (18) Gr as g$1(~T*/R)R Unfortunately, the details of the radial temperature gradient are extremely difficult to analyze without full solution of the field equations, especially when the boundary conditions allow for the

of the basic behavior of the flow and solute segregation, as noted by several other studies of convection in crystal growth [26—28].Two flow regimes arise most frequently in convection driven by lateral temperature gradients. We consider the limits for which thermal convection is either so weak that it can be analyzed as a perturbation from the uniaxial motion for crystal growth or motion so strong that the flow field has boundary layers. For small driving forces, i.e. small values of Gr, the about flow can be analyzed series written the static state of as Gr a= Taylor 0, as done in refs. [29,30] and by others for simpler thermal boundary conditions. For low Prandtl number fluids the temperature field is set entirely by conduction, and only viscous and buoyancy forces

162

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification

arise in the momentum balance at the leading order in the expansion. Then the intensity of the flow in this limit should be O(Gr) and is written in terms of the circulation ‘I’~as ~*(r* z*) = O(Gr). (20) The slow flows that occur when Gr ~u 1 lead to solute segregation only when they are strong enough to distort the uniaxial streamlines that correspond to diffusion-controlled crystal growth. The dependence of solute convection in this regime depends on the magnitude of the solutal Peclet number Pe~as Sc~Gr, which may still be extremely large because of the small mass diffusivity associated with solutes in the melt, and the magnitude of the growth velocity as set by the mass transfer Peclet number corresponding to the growth velocity PeSc/Pr. When momentum convection is weak (Gr ~ 1), but solute convection is strong (Pe 8>> PeSc/Pr> 1), the variation of the concentration field along the melt/crystal interface should scale as O(P~1/3) for a closed cavity flow near a no-slip boundary, as discussed by Pan and Acrivos [31]: 3= C(r) (Sc Gr) —1/3 c(r, h(r)) = C(r) pe~’/ (21) .

,

where C(r) = 0(1) for 0 ~ r ~ A. The dependence of the radial segregation on the intensity of convection should follow eq. (19). If a diffusion layer thickness 8 is defined by equating the total mass flux to an average concentration difference divided by this length scale, the dimensionless (scaled with the length of the ampoule) boundary-layer length will depend on O(Pe~1/3)~ The dependence of the effective segregation coefficient on 6 is derived in our notation from ref. [9] as keff = k + (1

—

k k) exp( a6 Pe Sc/Pr)’ —

(22)

As P; increases, 6 0 and keff —* k as O(Pe~1/3) or O(Gr1~3). The dependence of the flow on Gr becomes significantly more complicated when the Grashof —~

number is large enough that inertial effects and energy convection become important. Inertial boundary layers develop along solid surfaces, but interact across the cavity of melt and are not easily analyzed. Scaling for the dependence of the flow intensity on Grashof number is possible in the boundary layer limit where a nearly isothermal vertical boundary is isolated from viscous flow occurring along other surfaces in the cavity. Then the flow intensity ‘I” varies with Grashof number according to

~h1*(r*,z*)

=

O(Gr”2),

(23)

for low Prandtl number melts [32]. The variation of the concentration field along this vertical boundary then scales as -

c(A, z)

=

O(Gr Sc/Pr)

1/4

,

(25)

-

and the boundary layer thickness 8 scales as O((GrSc/Pr) ~-1/4) for Pe8>> PeSc/Pr and P;>> 1 [33]. Unfortunately, the flow near the melt/crystal interface in the vertical Bridgman geometry is not at all likeflow the along flow along a vertical butcavity. is the turning the bottom of awall, closed Only a few attempts [34] have been made to determine the proper scalings for flow boundary layers in such complex geometries. The closed streamlines for the bulk convection in the PSSM cause the value of the concentration in the uniform core to be a function of flow intensity. Since most solute is contained in the core, this value sets the effective segregation coefficient for intense convection; see eq. (1). Dimensional analysis of the concentration in the core predicts that it depends on RaT according to ccore = CO(RaT)~, where nk and C 0 are constants. 2’~) in From eq. (1), the effective segregation coefficient should vary with RaT according to O(Ra~’ the PSSM. The results for the low Grashof number and boundary-layer limits described above are only indications of the regimes for momentum and solute transport in the system studied below. The results described in sections 5 and 6 are compared to the scaling analyses reported in this section.

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification

163

2.3. Thermal boundary conditions for furnaces

3. Petrov—Galerkin finite-element analysis

The thermal boundary conditions for the MIT heat-pipe furnace are described by particular functions for the Biot number and ambient ternperature profile in eq. (15). The changes in the heat transfer coefficient between the three zones of the furnace are modelled by the function

The Petrov—Galerkin finite element formulation and the Newton iteration used for solving the steady state, free-boundary problem given by eqs. (1)—(17) are described fully in ref. [20]. For an approximate interface shape, the regions of melt, crystal and ampoule are divided into quadrilateral elements and Galerkin finite element approxima-

Bi( z)

where the constants in this expression have been fit to the measured wall temperature of the furnace of Wang and Witt [14]. The constant gradient furnace built by the Grenoble group is modelled by using a constant Biot number Bi = 100 over the entire length of the ampoule and specifying the furnace temperature profile in dimensionless form as

lions are constructed for the field variables. The temperature is approximated in all three phases by an expansion of Lagrangian biquadratic basis function. The velocity and solute concentration in the melt also are expressed in expansions of Lagrangian biquadratic functions. The pressure field in the melt is represented in terms of continuous bilinear Lagrangian basis functions, as suggested in the original mixed finite element method for the Boussinesq equations [36]. The interface shape h(r) is approximated by a series of Lagrangian quadratic polynomials and the coefficients in this expansion are computed as part of the solution of the entire free-boundary problem. This approach is an extension of the methodology put forth in ref. [37], which has proven extremely successful for solving solidification problems. Weighted residual equations for the field equations (6)—(11) are constructed using the Petrov—Galerkin formulation introduced by Zienkiewicz and Heinrich [38] for solving convectiondominated partial differential equations with finite

9 Iz~= 1

element approximations. weighting function used to form In this the residual methodequathe

=

C4 (1 + 0.5 [c1tanh C2 (z + tanh( c (z + c

))] }

—

2

3

—

C3) (26)

‘

as designed by Bourret et al. [17]. The coefficients (C1 have the following significance: C1 sets the ratio between the Biot numbers in the hot and cold zone, i.e. c1 as Bih/Bi,.; C2 sets the sharpness of the transition in Bi(z) between the isothermal zones and the adiabatic region; C3 is the dimensionless length of the adiabatic zone; C4 determines the absolute values of the Biot numbers in the hot and cold zones. The furnace wall temperature is modelled by the function 0 (z) = ~[i + tanh(12(0.5 z))], (27) }

—

—

z

(28~ /

The thermophysical properties of ampoule materials used in the calculations are listed in table 3 and are assumed to be constant over the temperature range of interest. Wang [35] argues that this is a good assumption for graphite and boron nitride. We consider only conduction through the ampoule wall and convection caused by its translation. Neglecting internal radiation through the wall is valid for boron nitride and graphite ampoules, but may be a serious error for systems using quartz. Perfect thermal contact is assumed between the ampoule and the crystal, so that the contraction of the solid during cooling of the ampoule is not taken into account.

tions depends on the magnitude and direction of the velocity field within each element, so that the directionality of convective transport is accounted for in the algebraic equation set. Petrov—Galerkin methods reduce to classical upwinding techniques used in finite difference formulations when applied to the simplest finite elements and to the convective diffusion equation with constant fluid velocity. The Petrov—Galerkin formulation used here is more sophisticated in that it reduces to Galerkin’s method in the limit of low elemental Peclet number, Pe* as V*he/a*, i.e. when either the flow speed (V *) or the element size (he) is small compared with the magnitude of the ap-

164

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification

propriate diffusivity a*. Two distinct velocity scales are present for solute transport: the transla tion rate measured by the Peclet number for mass transport scaled with the growth velocity PeSc/Pr and the velocity of natural convection which scales as (RaSc/Pr)” where the value of the exponent p is discussed in section 2.2 and Ra is either the thermal or solutal Rayleigh number. We have found that more accurate solutions for the solute field are computed when S(r, z) is decomposed into parts corresponding to each of these mechanisms; the details of this decoupling are given in [20]. We apply the Petrov—Galerkin method to the momentum, energy, and species transport equations, although the Prandtl number is small enough that no effect occurs from the extra weighting function in the energy equation for the convection levels computed here. The extensions of the methodology introduced in ref. [38] that are necessary for implementing the Petrov—Galerkin weighting functions with non-rectangular elements are discussed in ref. [20]. Each of the energy and species flux boundary conditions (eqs. (13), (14) and (16)—(18)) are incorporated into the residual equations for the field variables as weak boundary conditions. A separate set of residual equations for locating the shape of the melt/crystal interface is constructed from the melting point condition which is distinguished as the constraint for determining the interface, as suggested in ref. [37]. The entire set of nonlinear algebraic equations for the field variables and interface shape is solved simultaneously by Newton’s method as described in ref. [10]. Quadratic convergence is achieved when the dependence of the Petrov—Galerkin weighting functions on the velocity field is ineluded in the formulation of the Jacobian matrix at each iteration. Besides the rapid convergence rate Newton’s method is extremely useful as the basis for computer-aided perturbation methods for detecting multiple solutions to the equation set and for fully implicit time-integration methods for computing transients including the moving boundary; see refs. [39] and [40] for details. A typical finite element mesh used in the calculations for the MIT system is shown in fig. 2. The discretization used for most of the calculations contained 20 radial elements in melt and crystal

__________ -__________

_______ _______ ______

MELT

AMPOULE

MELT/CRYSTAL INTERFACE

CRYSTAL

I

it,, Fig. 2. Sample finite element mesh used for calculations in the MIT system.

and 4 elements across the ampoule. The axial approximation was composed of 48 elements in the melt and 16 elements in the solid. This mesh resulted in 26,600 algebraic equations. The Jacobian matrix is factored using a modified frontal algorithm [41] and its factors are stored using buffered, asynchronous input/output designed to minimize the number and increase the speed of data transfers. Each Newton iteration required 400 s on a Cray-iS supercomputer. Convergence

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification

was achieved in fewer than 6 iterations. Checks on the accuracy of the finite element results were performed using coarser meshes and are discussed in ref. [20] for an idealized growth system. Because of the high Schmidt number associated with the materials studied here the largest error occurred in the solute field. We report only calculations that satisfy the species mass balance at the melt/solid interface 1 ç —

,.

A

1

.,

c~r,h~r))rdras
i

=—

m

k

~29)

to within 10%. The calculations are reported in the form of contour plots of the stream function, temperature, solute concentration, and melt density throughout each region. The stream function is calculated by the scheme described in ref. [9]. Segregation of dopant across the crystal is reported in terms of the percentage of radial segregation defined as =

{

(

max [c r, h (r))] 0~r~A

—

0 ~mm r~A

[c(r, h(r))1}k.

(30)

Axial mixing of dopants is measured by the effective segregation coefficient defined by eq. (1) where the average are computed from the finite element solutions of the PSSM.

4. Effect of ampoule on idealized bridgman growth

system In the previous calculations for a vertical Bridgman system [9—12,20],perfect communication has been assumed between the furnace wall and the melt and crystal. This corresponds to assuming Bi —* oo and an ampoule thickness ta = 0 in the model presented here. Then the discontinuity in thermal boundary conditions at the junction of the adiabatic and hot zones of the Bridgman system leads to large radial temperature gradients near this point. The flows computed in each of these reports were driven primarily by this mecha-

165

nism. We have tested the effects of introducing an ampoule and a finite heat transfer rate between the ampoule and furnace on the magnitude of the axial and radial temperature gradients for the idealized system used in ref. [9]. These calculations are based on the set of thermophysical properties listed for gallium-doped germanium discussed in section 5 with the exception that the thermal conductivities of melt and crystal were assumed to be equal at the value appropriate for the melt. The ampoule is assumed to be four times longer than .

.

its radius, i.e. A = 0.25, the ratio of the thermal conductivity of the ampoule to the melt and crystal is taken as Ka = 0.5 when it is included. The radial and axial temperature gradients predicted for pure conduction (RaT = 0) are shown in fig. 3 for (a) the idealized systems used in ref. [9], (b) the same system with a large heat transfer coefficient, Bi = 100, and an ampoule of thickness Ia = 0.125, and (c) a realistic set of heat transfer coefficients (Bi = 1.0 in hot zone, 0 in adiabatic zone, and 0.5 in cold zone) and an ampoule. The radial temperature difference is defined as 9R = [0(z, 0)— 0(z, A)]/9(z, 0). (31) ~

The large non-uniform radial gradients in the idealized system are apparent. Introducing either finite heat transfer coefficients, the thickness of the ampoule, or both of these effects decreased these gradients along with the magnitude of the axial gradient. When both the finite heat transfer coefficient and the ampoule are introduced the radial difference ~°R dropped from 0.31 totemperature 0.08. The intensity of the convection decreases as the magnitude of the radial temperature difference is diminished. This effect is shown in fig. 4 by the value of sI,* associated with the convection cell in the idealized system, with and without the added thermal resistances caused by the ampoule and the finite heat transfer rates between the ampoule and the furnace. At RaT = 1.0 X iO~,the intensity of the cellular flow has decreased by more than an order of magnitude. As a result, quantitative prediction of the flow field and dopant segregation requires accurate calculation of the heat transfer in the furnace—ampoule system.

166

P.M. Adornato, R.A. Brown

-0.35 0.9

0.8

—

/ Convection and segregation

DIMENSIONLESS RADIAL TEMPERATURE 8R DIFFERENCE ~ -O 5 - 0.05 0.05 0.15 0.25

-0.25 I

I

I

I

in directional solidification DIMENSIONLESS RADIAL TEMPERATURE DIFFERENCE AØR

0.35

-025

-0.15

-0.05

0.05

0 15

025

b

—

°~

MELT No Ampoule B~--m

MELT Ampoule 08 ‘S

La

.5’

40.7-

La I-.

-

-

sfO7

55.

z a

—

‘.5

2 -

‘5’,

o

a

)

0 006_

-

/

—— -J 4

~05

-

40.5

-

/ La

2 (I) z La

0.4

-

r’ ‘5’ •5.

0.4

-

—

0 — (0

z

•.. ‘S.

a

~0.3-

-

\

a

\ \

a

0.2—

0.2-

-

CRYSTAL 0.1

I

I

I~

I

CRYSTAL

I

0.1 I I I 0.0 02 04 06 0.8 10 DIMENSIONLESS CENTERLINE TEMPERATURE

0.0 0.2 0.4 0.6 0.8 1.0 DIMENSIONLESS CENTERLINE TEMPERATURE

DIMENSIONLESS RADIAL TEMPERATURE DIFFERENCE -0.15 -0.05 0.05 0,15 0.25

-0.25 05,

08

—

I

MELT Ampoule Bi’OlI)

I

I -

Lu 2

a -

//

-J 4

0.5

-

// I

(0 U)

La

0

I

z La

I

a

-

C

0.2—

-

CRYSTAL

I

0.I I 00 0.2 0.4 0.6 0.8 1.0 DIMENSIONLESS CENTERLINE TEMPERATURE

Fig. 3. Axial temperature profile and radial temperature differences for vertical Bridgman furnaces with varying thermal boundary conditions and ampoule configurations. The three cases correspond to: (a) system without ampoule and large Biot number; (b) system with ampoule and large Biot number: (c) system with ampoule and Bi = 1.0. The axial coordinate is measured from the bottom of the ampoule.

P.M. Adornato, R.A. Brown

I

/ Convection and segregation

-

167

typical of experiments by both the MIT and Grenoble groups do not lead to solutal-driven

I 0 WITHOUT THERMAL RESISTANCE O WITH THERMAL RESISTANCES

.01

in directional solidification

.

convection. This conclusion is based on the magnitude of the solutal Rayleigh number Ra~= 0(102) estimated to be appropriate for this concentration and the results for the thermosolutal

._—“~

.001

0001

L 102

I

I

IO~ THERMAL

Io~

Io~

RAYLEIGH NUMBER

Fig. 4. Effect of heat transfer resistances on the intensity of the cellular convection without and with the added thermal resistances of the ampoule and finite heat transfer coefficients between the furnace and ampoule for the idealized growth system used by Chang and Brown [9].

5. Growth of the gallium-doped germanium The thermophysical properties used in this study for gallium-doped germanium were compiled from results of Crouch et al. [45] and Wang [35] and are listed in table 1. All properties were assumed to be independent of temperature. This assumption is substantiated at least in part by the data in ref. [44]. We assume that the small concentrations of gallium (O(10i8) atoms/cm3 or O(10~) mol%)

The calculations for GaGe growth are presented separately for the two furnaces and the different ampoule designs. The thermal Rayleigh numbers appropriate for solidification on earth

are reported in table 2. The Petrov—Galerkin/ finite-element method easily computed the temperature and flow fields for the growth of dilute germanium under these conditions. Temperature and velocity fields were computed for RaT = 3 X 108, the value appropriate for the MIT furnace. However, the resolution of the mesh was insufficient to resolve accurately the boundary layers present in the solute for the most intense flows. Accurate solute fields were only computed with this mesh for RaT < 5 X iO~. 5.1.

Vertical-Bridgman heat-pipe furnace.’ boron

nitride ampoule 5.1.1. Temperature and flow fields Calculations for growth of GaGe alloy in a

Table 1 Thermophysical properties of germanium Thermal conductivity of the melt, Km (W/°Ccm) Thermal conductivity of the solid, KL (W/°C cm) Density of the melt, Pm (g/cm3) Density of the solid, p~(g/cm3) Specific heat of the melt, Cpm (J/°C.g) Specific heat of the solid, ~ (J/°C g) Melting temperature, Tm (°C) Slope of the liquidus curve for dilute silicon, m (~ C/mol% . Si) Kinematic viscosity, v (cm2/s) Heat of solidification, LI !I~(J/g) Thermal expansion coefficient, f~(°C~) Solutal expansion coefficient, $ (mol% Si)~’ Mass diffusivity of Ga in Ge, D (cm2/s) Equilibrium segregation coefficient of Ga, k Equilibrium segregation coefficient of Si, k

0.39 0.17 5.5 5.5 0.39 0.39 937.4 11.3

o.ooi3 460

5.oxio4 —

2.1 x i0~ 0.087 4.2

168

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification

Table2

R

Dimensionless groups and their charactenstic values for mod elling gallium doped germanium and sihcon—germamum crystal growth

I

0

u

0

3

________________________________

Thermal Rayleigh number

RaTm g$1 LITL

Solutal Rayleigh number

Ram

g$ c

-~ I-

VL Pe~~t__

Prandtl number

PrU)

_____

Oto—10X107 MELT

16x10’2

:5:5,1

AMPOULE

________

—

Sc

a -~

70x103

MELT

62 INTERFACE

St

pç, LIT

“

N

11

LI/I Stefan number

*55

Oto3OXlO

3 0L

Dimensionless growth velocity

Schmidt number

ROT 5

T ii

~

INTERFACE

~

40 CRYSTA

boron nitnde ampoule were performed for several combinations of growth speed gradient zone length and ampoule thickness for the hot and cold zone temperatures most used in the experi mental system Parameters corresponding to the base case for this system are tabulated in table

CRYSTAL

——

3.

‘~‘

Heat transfer in the heat-pipe system differed dramatically from the idealized system described in ref. [9] and in section 4 above. The combination of the long gradient zone and the different, but comparable, thermal conductivities between melt, crystal and ampoule led to a region of melt adjac-

_________

4 Fig. 5. Sample temperature fields for growth of GaGe in MIT furnace with thermal Rayleigh numbers of (a) 1 x i05 and (b) 5 x iO~Ra, = 0. Results are for boron nitride ampoule.

Table 3 MIT and Grenoble Furnace Design Parameters Parameter Ampoule length, L (cm) Crystal radius, R (cm) Ampoule outer radius, Ra (cm)

MIT ampoules

Grenoble ampoules

Boron nitride

Quartz

Graphite

Composite

7.62 0.63 0.95

7.62 0.83 0.95

7.0 0.50 0.70

7.0 0.35 0.50 (0.70)

3.81 300 4.0

3.81 300 4.0

—

—

350 0.83

350 4.0

0.26

0.028

3.26

3.26

Gradient zone length, L 5 (cm) Temperature difference, T~,— 7’~(°C) Growth velocity, i~(gem/s) Thermal conductivity of ampoule, ka (W/°C~~m)

(0.028)

P.M. Adornato, R.A. Brown

.25

I I

I

/

ent to the interface where the fluid near the centerline was hotter than material near the ampoule. An interface shape concave to the solid was always computed regardless of the position of I

ROT’ I a IO~

Rag=0

~ uS

2

-

v9

PIPE GRENOBLE HEAT FURNACE

the ampoule inside the gradient zone. Temperature fields for RaT of 1 X iO~ and 5 X iO~are presented in fig. 5 and show these characteristics. The large convection level associated with the higher value of RaT did not deflect the isotherms

C ONSTA~TN~ADIENT I .2 .4 .6 .8 DIMENSIONLESS RADIUS r

Fig. Bridgman 6.= 1Interface (MIT)The shapes andradial constant computed for GaGe (Grenoble) growth have infurnaces; vertical RaT xtoiO~. and gradient axial coordinates resealed account for the differences in the inner radii ofbeen the two ampoules.

‘4~m/s

.15

.1

-

.050

169

Convection and segregation in directional solidification

Streamlines 6

~ STREAMLINES

r

Ra1’lxlO

RATS

lx5 IO

Streamlines Ra

—

1 ‘5xI0~ R~

—

7 1’IxIO

I-

o

0

I

I

*~

I-

i~,

2.93 x I0~

4

MAX

5’987e I0-~

0 I

2

I-

‘6.22xI0

‘~‘mox

/ /~

6.22s IO~

I..

8x IO~

C

IO3

/ S C nix’

° N

—6.7

‘0

Ui

0

0 Z

‘4

0

~ Interface 014

a

_ _

INTERF~E

u4mls’ -2.26 — 4

N

x l0~

m

(~~‘ s

t~FACE _

/“

4 a — —1.59

// ‘~

‘a

IO~

1~

~ _*,~lnterface 0 5,5~_4.19~t0A

2

Fig. 7. Sample flow fields for growth of GaGe in MIT system with boron nitride ampoule. Thermal Rayleigh numbers are between 1 x 106 and 5 x 10~.Streamlines are spaced at equal intervals between the maximum (or minimum) values for the cells and zero.

170

P.M. A dornato, R.A. Brown

/

Convection and segregation in directional solidification

appreciably. The corresponding interface deflection is due mainly to the differences in thermophysical properties and to latent heat release. The interface shape for RaT = 1 x iO~is plotted as fig. 6. The changes in the thermal field caused by the ampoule led to markedly different melt flows from those predicted in the idealized system. Instead of the single toroidal cell flowing down along the centerline and up along the ampoule wall, there were two axially stacked counter-rotating cells, with the cell nearest to the interface moving downward along the wall and up at the centerline. Representative contours of the stream function for these flows as a function of RaT are shown in fig. 7. For low values of RaT, the intensity of the upper cell is greatest; however, the intensity of the cells becomes almost equal as the Rayleigh number is increased. This result is in contrast to those found for the prototype system [9,20] where secondary cells formed only as the result of intense circulation of the primary cell and, hence, were much weaker than the primary flow. The intensity as measured by the absolute values of the stream function ‘I” for each flow cell is plotted separately in fig. 8. The plateau in ‘I” for the upper cell at low values of RaT represents the volumetric flow rate for the uniform motion caused by the crystal growth rate J’~.The increase of the flow intensity with Rayleigh number is approximately linear up to RaT = 1 X iO~,indicating that inertial forces are not important and the scaling described by eq. (20) is valid. The top cell seems to develop a boundary layer structure for RaT> 1 x iO~,indicative of the slower increase in the circulation rate with RaT. The slope of the curve should become roughly one half if the boundary layer scaling for a vertical wall, eq. (23), is applicable. The value of the slope at the highest value of RaT agrees approximately with this prediction; calculations for more intense flows are needed to make the comparison precise. Although of comparable intensity, the convection in the lower cell has not developed boundary layers for the same values of RaT. This is expected since the lower cell represents a turning flow near the melt/crystal interface and is much more restricted by the geometry than the top cell.

102

I

R O~ V 9

I

I

I

I

o

‘4j~mis

IO-~ -

-

UPPER FLOW

10

—

-

LOWER FLOW CELL

I O~ -

-

I I I I IO~ IO~ 106 IO~ 108 THERMAL RAYLEIGH NUMBER ROT Fig. 8. Maximum value of stream function in the top and bottom toroidal cells as a function of RaT for growth of GaGe in MIT system with boron nitride and quartz ampoules. I ~6

IO~

5.2. Solute segregation 5.2.1. Solute fields

The two-cell structure of the flow in the vertical Bridgman furnace has a profound effect on the structure of solute segregation in the melt and crystal. Solute profiles for dilute gallium in germanium are shown in fig. 9 as calculated for the flows shown in fig. 7 and a crystal growth rate of 4 tim/s. The solute field is approximately uniform radially for RaT ~ 1 X iO~.Increasing the level of convection sweeps solute near the interface from the wall of the ampoule toward the centerline of the melt and leads to higher dopant levels near the center of the crystal. The motion of the upper cell sweeps solute in the opposite direction. As convection is increased the melt within the core of each cell becomes better mixed and gradi-

P.M. A dornato, R.A. Brown

[~

Ac

•IxIO~ 1.80

I I I

II

/

lAc

I I

0

I

I-.

171

Convection and segregation in directional solidification

1.92

I I I

I

I R 0T

T

uIxIO Ac ‘1.197

I

I0

I

I

I

I

I I

I I I

I0

x

I—

4.86

I

I

1w

I I

I I

12.84

I~

17.24

,______/“

l~

L__—=—1 N

3.38

g

0 N

N

F—=.....

4.

21.02

I

.~

4. 12.95 INTERFACE

Fig. 9. Sample gallium concentration fields for growth of GaGe in MIT system with boron nitride ampoule; 0 ~ RaT ~ 1 x ~ = 4 pm/s.

ents of solute concentration are confined to the

boundary layers that define the edges of the cells. These layers are not uniform at RaT = 1 X iO~, leading to substantial radial segregation in the crystal. Moreover, the region of weak convection that separates the two cells supports an internal layer in the solute field. The formation of the internal layer complicates the numerical solution of the solute problem, because high accuracy of the finite element approximation is needed there as well as along the boundaries. The solute concentration in the melt, along the

and

interface is plotted in fig. 10 as a function of RaT. The radial segregation present for the low convection intensity (RaT = 1 X 10~)is caused mostly by the curvature of the melt/crystal interface. As in the prototype system [9,20], there is a maximum in the radial segregation of dopant which occurs in the boron nitride ampoule at roughly RaT 1 X 106. The radial concentration variation seems to approach a uniform value over the bulk of the interface with variations confined to the regions near the centerline and the ampoule wall. This profile shape is expected for intense convection

172

P.M. Adornato, R.A. Brown

/ Convection and segregation I

I5—.~~

in directional solidification

I

I k ‘0.087

RaTrlxIO6

4— I

VgC4~.i.M/S

—

x IO~ —

~

1x105

I

~I2-

0

I

I

.020092

0

I

.040184

DIMENSIONLESS

I

.060276

.080368

.I0046

RADIAL COORDINATE r

Fig. 10. Variation of gallium concentration along the melt/crystal interface with RaT for growth of GaGe in MIT system using a boron nitride ampoule; ~ = 4 pm/s.

15

I!II~III~IIII~II

Ra

VERTICAL BRIDGMAN FURNACE GaGe

7 1:1xI0

cc

Vgr4~.m/s

‘\

-

\Ra~ \1x106

-

z

C I4 IZ5~r~

LU

o z 0

-

IxIO 5

0

0 0.0

I

I 0.2

I

I 0.4

I

II

I

0.6

I

I 0.8

1.0

NORMALIZED AXIAL COORDINATE, ~ Fig. 11. Axial profiles of gallium concentration for three values of RaT in MIT system with boron nitride ampoule.

P.M. Adornato, R.A. Brown

/ Convection and segregation

where the thickness of the solute boundary layer becomes vanishingly small everywhere except in the turning regions. The approximate boundary layer thicknesses predicted for these calculations is

V

9=2~.cm/s

discussed below in terms of the effective segregation coefficient. The effect of the multi-cellular convection on the solute field is most evident from axial profiles of the solute concentration along the centerline of the melt; these are shown in fig. 11 for three values of RaT and = 4 tim/s. The profile for RaT = 1 x iO~is essentially exponential. Two regions of almost constant composition separated by the steep internal layer form as the Rayleigh number is increased.

173

in directional solidification

AC ‘0 SI

V9’81.cm/s

AC~I 57

~

5.2.2. Effect of growth rate

The level of radial segregation depends on the interaction of the cellular convection with the growth velocity V~and the scale of the exponential concentration layer adjacent to the melt/ crystal interface. Calculations were performed at growth rates of 2, 4, and 8 tim/s to explore this effect. The solute fields for RaT = 1 x iO~and the lowest and highest of these growth rates are shown in fig. 12. The concentration variations in the lower cell are compressed somewhat by the thinner solute boundary layer at the interface associated with the higher growth rate at a fixed Rayleigh number. the solute layer inhibits mixing ofThe the thickness solute andofleads to larger variations of composition along the interface, as demonstrated by the interfacial solute concentrations shown in fig. 13. The amount of radial segregation defined by eq. (29) is plotted in fig. 14 as a function of RaT for the three growth rates. In each case, the maximum in the segregation occurs at an intermediate value of the Rayleigh number. The magnitude of the radial segregation scales almost linearly with the growth rate. Axial segregation of solute can be estimated from the results of the PSSM by computing the effective segregation coefficient defined by eq. (22). These values are reported in table 4 for the three growth rates. The decrease in keff toward k = 0.087 is obvious for all values of J’. The most effective axial mixing is found for the lowest growth rate, as suspected from the behavior of the radial segregation with increasing RaT.

7.11

___________

.~s

.81

0,55

11.98

5 27

Fig. 712.plots Sample gallium concentration fields for growth of are for growth rates of (a) V~ = 2 pm/s and (b) GaGe nitride ampoule and RaT 1 =8 pm/s. x i0 in MIT system with boron

5.3. Heat-pipe furnace.’ quartz ampoule 5.3.1. Temperature and flow fields

Using a quartz ampoule decreases the conductivity of the ampoule to a factor of 14 lower than the value for the melt and leads to the ampoule absorbing more of the radial temperature difference necessary to transport heat in the system. The temperature field computed for the quartz ampoule is shown as fig. iSa. The radial temperature gradients near the interface were lower, the interface was flatter (see fig. 6), and convection near the interface was weaker than for the boron

174

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification ‘‘‘I’’’

‘I’’

1I’’’1I

‘

III

1 111T11

U

VERTICAL BRIDGMAN FURNACE GaGe 7 ROT’ Ix10

0

Q

III

-

I‘ax

~ 10

-

I-

z

LU (_.)

-

z 0 U

U, 5U)

LU

Vg ‘4~m/s

-J

z

0

-

Vg’8~m/s

a LU

~ ~

0 0.0

I I

I

I

I

I I I I I I I I 0.2 0.4 0.6 0.8 NORMALIZED AXIAL COORDINATE,

e

I

20.5

I

I I I

I 1.0

I k ‘0.0877

18

—

RATC I x10 V 9’

5.5

-

8j.~

4MM/S 2

I

0.5

b 8 0

I I I I .020092 .040184 .060276 .080368 DIMENSIONLESS RADIAL COORDINATE r

.10046

Fig. 13. Variation of gallium concentration along (a) the centerline of the melt and (b) the melt/crystal interface with Vg for growth of GaGe in MIT system using a boron nitride ampoule; RaT = 1 x iO~.-The normalized axial coordinate in (a) is defined as

Em[z—h(O)]/h(O)

P.M. Adornato, R.A. Brown

/ Convection and segregation in

interface shape during growth in a quartz ampoule, while a concave interface is predicted from our

100

z

semi-transparent ampoule.

(.9 LU

~r 60 (9

0 80 /~N\\

calculations. This difference is probably a result of 5.3.2. Solute fields

.

LU

C,) —i

175

directional solidification

20 40

~

i0~

2

I05

RAYLEIGH

7thermal arewith shown in fig. our Solute neglecting fields computed for transport through Rayleigh the 16 ampoule. for theof vertical The concentration Bridgman system for the quartz lower value of RaT is uniform, that solute numbers 1 xradiative 106 andheat 1 xindicating i0field

- - -- -.-~.

06

IO~

NUMBER

08

Fig. 14. Percentage radial segregation as a function of thermal Rayleigh number for growth of GaGe in MIT system with boron nitride ampoule; results are shown for growth rates ~ of 2, 4, and 8 pm/s. Solid dots correspond to experimental measurements of Wang [35].

nitride ampoule discussed above. Large radial temperature gradients existed near the intersection of the adiabatic and hot zones. These gradients induced an upper flow cell that was more intense than the motion in the boron nitride ampoule. Flow fields for RaT equal to 1 x 106 and 1 X iO~ are shown as figs. i5b and i5c. The disparate strengths of the two cells is apparent here and from the values of the circulations within the cells which are plotted in fig. 8. The qualitative agreement between these calculations and experiments is not as close as for the results with the boron nitride ampoule, discussed in sections 5.1 and 7. Wang [35] observed a convex

transport is nearly diffusion-controlled. For RaT = 1 x iO~a core of constant concentration has not formed in the lower cell and the radial segregation for this case is higher than for RaT = 1 X 106, see fig. 17. Both results suggest that the radial segregation is governed by the incomplete mixing in the lower cell and not by the interaction of a boundary layer along the interface with the bulk flow. The high circulation rate in the upper cell causes axial segregation to be more severe in the quartz ampoule than in the boron nitride ampoule, even though the mixing in the lower cell is less effective. This is suggested by noticing that the convection of solute into the melt above the adiabatic zone is present in the quartz ampoule (fig. 16a), but is absent in the boron nitride system (fig. 9a). 5.4. Constant gradient furnace: graphite ampoule 5.4.1. Temperature and flow fields Typical design parameters used in the constant gradient furnace built in Grenoble are listed in

Table 4 Effective segregation coefficients computed from eq. (1)

MIT: Boron nitride ampoule

Grenoble: Graphite ampoule Grenoble: Quartz/graphite composite ampoule

RaT=1X105

RaT=1X106

RaT=1x107

(pm/a)

<(c))

~

ac))

~

ac))

keti

2 4 8

3.69 2.34 1.65

0.27 0.43 0.60

5.99 3.81 2.46

0.17 0.27 0.41

7.35 5.14 3.32

0.14 0.19 0.29

4

2.57

0.38

3.53

0.28

5.00

0.19

4

2.49

0.41

2.47

0.40

4.72

0.21

176

P.M. A dornato, R.A. Brown

/

Convection and segregation in directional solidification

Rat’ Ix 106 AMPOULE

5

i:

‘i:

MELT

~

X

I0~

I-.

0 I

I-

~I

INTERFACE

R~’ I

0 I I.-

I— 0

I I—

/ z 0

:.

—

~max

U

_______

~max

3

I

+I63x10

‘4’min —2.56x10

______

CRYSTAL

I

. mm

783x10

—

-4.48xI0

4.

~

~“INTERFACE

~-

~INTERFACE

Fig. 15. Sample temperature and flow fields for growth of GaGe in MIT furnace with the quartz ampoule; Ra~= 0. The plots are: (a) the temperature field for RaT = 5 x iO~(b) flow field for RaT —1 x 106; (c) flow field for RaT = 1 x iO~.Streamlines are spaced at equal intervals between the maximum (or minimum) values for the cells and zero.

table 3. A sample finite element mesh for the calculations used in analyzing this furnace geometry is shown as fig. 18a. When a graphite ampoule with high thermal conductivity is used the radial gradients are very low except near the melt/crystal interface, where large differences in temperature still exist because of the mismatch in the thermal conductivities of the three phases. The isotherms appear almost flat on the scale of the plot shown as fig. 18b, indicating the very low driving forces for convection in the system. The deflection of the

melt/crystal interface is greater in the constant gradient furnace, as indicated in fig. 6. Eliminating the large radial gradients that arise at the intersection of the hot and adiabatic regions in the heat-pipe furnace removes the driving force for the upper flow cell seen there; this difference is evident in the flow fields shown in figs. 18c and i8d where only the flow cell nearer the melt/ crystal interface has formed. As in the MIT furnace, this cell transports melt up along the axis of the ampoule and down along the wall. The inter-

P.M. Adornato, R.A. Brown

ROT’ 1xI06

RoT’ I~IO~

Ac-103

AcI.15

Lo 302

— ________

‘1 -0 -5I-.4 -at

__—~-‘~~~

/

IS ~

_________

ç

2.49

INTERFACE

Fig. 16. Gallium concentration fields for growth of GaGe in MIT system with quartz ampoule; (a) RaT = 1 x 106; (b) RaT =lx i0~.

I

I

Ro~0

____

2I3 Ra z~ 14 II910________________________________ RaT~IxI06 ~I0I -- V I I I 0.022 0.044 0.066 0.088 0.. I

4 IU 0 I0 — Ll -~.

I-. — 0 -~ -‘6 - 4

I

INTERFACE

177

Convection and segregation in directional solidification

DIMENSIONLESS RADIAL COORDINATE 9’4~m/s

Fig. 17. Variation of gallium concentration along the melt/crystal interface for growth of GaGe in MIT system using a quartz ampoule; V 5 = 4 pm/s and RaT = 1 x 106 and lx io~.

-~

08 -~

R018 I x iO~

7

3.86

X IO—~

ROT:Ix10

MELT

______

___

~‘mIn

____

4.44 x 10-6 (~ MELT

AMPOULE 02

4

~mIn = 10 —I.04x

4.

AMPOUL E CRYSTAL

t Fig. 18. Sample temperature and flow fields for growth of GaGe in Grenoble constant gradient furnace with a graphite ampoule; plots are: (a) the finite element mesh; (b) the temperature field for RaT = 1 X io~(c) flow field for RaT = 1 x io~(d) flow field for RaT = lx i07.

178

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification

RaT= I x 106 Ra1r I x IO~ Ro1=Ix l0~ 2.34

4 75

2.I94

4.562— 8 II4 -~ 11.666—

8.89

~

4.

84I_

~

L

-

12.82-

12.08

4. 5 ~ RaT ~ 1 x l0~.

Fig. 19. Sample gallium concentration fields for growth of GaGe in Grenoble furnace with graphite ampoule;

face is flatter and the flow cell weaker for a given value of RaT than for growth in the MIT furnace with a boron nitride ampoule. The intensity of the circulation, ~I,*,scaled approximately linearly with increasing RaT for values up to 1 X iO~(see figs. 18c and 18d), the largest value reported. This is also true for the lower cell in the MIT furnace, as shown in fig. 8. Also, the rounded character of the

ix i0

gradients extend from the melt/crystal interface into the bulk melt and are axially aligned in the nearly uniaxial melt flow created in this region by ampoule translation. The radial segregation predicted for a given 50 I

~

recirculation for the highest value of the thermal Rayleigh number indicates that the dominant driv-

I

I

Ro

40

-

5 ‘0 Vg’4/.Lm/s

VERTICAL BRIDGMAN FURNACE BORON NITRIDE AMPOULE

4

ing force for the motion is the radial gradients near the interface, not non-uniform gradients along its length.

30 U ‘I, -J

5.4.2. Solute fields The concentration fields for gallium in the melt are shown in fig. 19 and correspond to the flow fields shown in fig. 18. Convection of the gallium becomes noticeable at RaT = 1 x 106; the maximum value of the interfacial concentration of galhum occurs along the centerline, as was the case for the MIT furnace. A uniformly mixed core begins to appear for RaT = 1 X iO~,but is confined to the region near the melt/crystal interface where convection is intense. Appreciable solute

4

20 ~APHITEAMP0ULE CONSTANT GRADIENT FU

I— U

0

w I

CONSTANT GRADIENT FURNACE: COMPOSITE AMPOULE I I0~ I 08

THERMAL RAYLEIGH NUMBER Fig. 20. Percentage radial segregation as a function of RaT for MIT (boron nitride ampoule), Grenoble (graphite ampoule) and Grenoble (composite ampoule) systems.

P.M. Adornato, R.A. Brown

/

179

Convection and segregation in directional solidification

1.0

09

IS

06

-

MELT 07 QUARTZ LINER GRAPHITE AMPOULE

06

~‘,I ‘ ,

GRAPHITE AMPOULE QUARTZ LINER

05 INTERFACE 0.4

c:’

CRYSTAL

02

01

II:’

0

RaT Ra

1

Ix

=

I x 106

I0~

RaT= lx

tkmax

2.25 x I0~

9

tI~lmaxxIQ5

—5.96 ~

io9~fl.$ 4..

—9.83

x

I05

Fig. 21. Sample (a) finite element mesh, (b) temperature field, and (c)—(e) flow fields for growth of GaGe in Grenoble furnace with composite ampoule.

180

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification

value of RaT is plotted along with the data for the MIT furnace in fig. 20. The value of & is comparable for both systems, as expected since the circulations near the interface are similar in form and magnitude. Interestingly, the axial segregation for the vertical Bridgman and the constant gradient furnaces also are predicted to be similar (compare values listed in table 4 for J’ = 4 rim/s) even though the axial mixing supplied by the upper cell in the MIT system is absent in the constant gradient system.

ampoule made of a quartz sleeve enclosed in a graphite liner; the dimensions of the ampoule parts are given in table 3. A typical finite element mesh, temperature field and flow structures for this ampoule are shown in fig. 21. Introducing the quartz sleeve reduced the radial temperature differences, near the interface and flattened it, as compared to the results for the graphite ampoule in the Grenoble furnace and the boron nitride ampoule in the MIT furnace. The intensity of the convection decreased by almost an order-of-magnitude from the value for the graphite ampoule; compare figs. 21e and 18d. In addition, the recirculating cell in the composite ampoule was not limited to the melt,~near the melt/crystal interface, but spread upward for increasing values of RaT. The flow in the lower portion of the ampoule has the form of a nearly one-dimensional motion coupled to the turning region at the interface. The flow pattern away from the interface resembles that of the parallel flow between the differently heated plates where the convection is driven en-

5.5. Constant gradient furnace: composite quartz/ graphite ampoule

5.5.1. Temperature and flow fields The radial gradients near the melt/crystal interface can be lowered further by introducing a more insulating ampoule material, as demonstrated for the case of quartz in the MIT system. Following the experiments of Rouzaud et al. [15], we performed calculations for a composite

Ra7’ I

x I0~

Ra7=

I

x 106

RaT

I

x

l0~

5.27

9.4e 4.21

4.26 7.52 10.78

I0.62

4-

4_

4.

Fig. 22. Sample gallium concentration fields for growth of dilute GaGe 7. in Grenoble furnace with composite ampoule; ~ ix i0

tx l0~~

RaT

P.M. Adornato, R.A. Brown

/

because of the stabilizing effect on the density associated with the silicon concentration in the

_____________________________________

Ra

I

I

I

I

I

V~’ 4~m/s

is

z

4

-

3

-

2

-

I I

-

I

-

~ 10 I— 0

‘I a I o~

Ra Ra ‘Ix i06

W 2

~

T

T

.

5

RaT’ Ix io I

I

I

181

Convection and segregation in directional solidification

I

0.01 0.02 0,03 0.04 0.06 DIMENSIONLESS RADIAL COORDINATE r

Fig. 23. Variation of gallium concentration along the melt/crystal interface as a function of RaT for growth GaGe in Grenoble system with the composite ampoule; V5 = 4 pm/s.

tirely by the resulting lateral temperature gradient [43]. 5.5.2. Solute fields Gallium concentration fields that correspond to the convective patterns given in fig. 21 are shown in fig. 22. Convective mixing is noticeable in the results of the two higher thermal Rayleigh numbers, but no well-mixed core is formed. The solute concentration along the melt/crystal interface for these cases is plotted in fig. 23. The radial segregation & is plotted as a function of RaT in fig. 20, where it is compared directly with the results for the graphite ampoule; & has reached a maximum in the graphite ampoule near RaT = 1 X iO~,but is increasing with increasing RaT in the composite ampoule. The effective segregation coefficients computed for the composite ampoule are listed in table 4 and are higher than the values for the graphite ampoule. This is indicative of the weaker convection in the composite ampoule.

6. Growth of silicon—germanium alloy: constant gradient furnace Introducing silicon as a solute into germanium has been suggested by Rouzaud et al. [15] as a mechanism of reducing convection in the melt

diffusion layer at the interface. Experiments in the Grenoble furnace substantiated this idea by showing the tendency of the segregation to evolve toward diffusion-controlled growth with increased silicon concentration. Introducing the silicon solute into the melt has other interesting implications for the flow pattern. The diffusion layer establishes an axially stabihizing concentration gradient that decreases with distance away from the interface while flow is driven by the radial temperature gradients along the axis. . . We expect that the sideways diffustve instability described by Hart [42] is possible. Hart analyzed the transition of a parallel flow between two vertical flat plates that establish a lateral temperature gradient which drives a flow in the presence of a constant axially stabilizing concentration gradient. He showed that the flow pattern undergoes a transition to a pattern of vertically-layered cellular flows as the solutal Rayleigh number is increased at constant RaT. When the mixing in these cells is intense the concentration field becomes nearly constant in each cell with internal solute layers separating the core values in adjacent cells. These cells have been observed by many researchers; see ref. [44] for references. We expect that flow patterns similar to those described by Hart are possible in the Grenoble SiGe experiment, but that the non-uniformity of the radial temperature gradient along the ampoule, the exponential shape of the axial concentration profile, and the ends of the ampoule all break the symmetries inherent in the idealized problem he studied. These imperfections cause the flow transition to the vertically-stacked cellular flows to occur smoothly with increasing solutal Rayleigh number Ra and force the cells to be localized in the concentration gradient near the interface. We have performed calculations for the growth of SiGe alloys in the graphite ampoule described in section 5.4 to test this hypothesis. We model the SiGe alloys as having the same thermophysical properties as germanium except for changing the segregation coefficient, liquidus slope, and diffusion coefficient to the values appropriate for silicon; these properties are listed in table 1. The first .

182

P.M. A dornato, R.A. Brown

/ Convection and segregation

set of calculations is for the value of the liquidus slope m = 0, so that the phase diagram does not affect the shape of the melt/crystal interface. The temperature fields for these calculations are almost identical to the case shown as fig. 20a, because the thermal properties have not been altered by the addition of silicon. Varying the solutal Rayleigh number Ra5, while holding RaT constant, is equivalent to increasing the concentration of silicon in the melt for a system with all other parameters held constant. When the concentration of silicon is low enough there is no interaction between it and the density of the melt, i.e. Ra = 0, and the velocity field will be identical with that predicted for the dilute GaGe mixture in section 5.4 for the same thermal Rayleigh number. Concentration fields for a dilute silicon solute corresponding to these cases are shown in fig. 24 and show the same trends discussed for the galhum solute, but with the concentration lowest near the interface because silicon’s segregation coefficient is greater than unity, For RaT ~ 1 x iO~,introducing the stabilizing

in directional solidification

solute field adjacent to the interface by decreasing the solutal Rayleigh number Ra from zero, diminishes the recirculation caused by the mismatch in thermal conductivities until the flow field is essentially the uniform axial motion corresponding to steady crystal growth. Increasing RaT to 1 x 106 intensifies the recirculating flow in the absence of the solute. The effect of adding silicon on the flow is demonstrated in fig. 25. Low concentrations (Ra5 ~ —5 X 10~)weaken the recirculating flow, but do not change its structure. The flow field for Ra5 = 1 x 106 shows undulations just above the recirculating cell driven by thermal convection which are indicative of the beginning of new flow cells formed by the sideways diffusive mechanism of Hart [42].These new axially stacked cells are more apparent in the results for Ra = 1 x iO~and have roughly the same intensity as the flow induced by crystal growth. Attempts to compute the solutions for values of Ra beyond 1 x iO~all failed indicating that no solution of the discretized equation set is possible. Since we have not attempted to use finer finite element —

—

—

0.5939 0.6662 _____

RaT’ I x I0~

6

0.6426

Ra

Ru T =IxI0 Ru 50

0.4690

______

Ras’0

T

1’IxIO Ra5~0

0.5007 05199 0.3440

_______

Q3589

0.3736

0.2274

4-

0.2190

02170

4-

t~.

Fig. 24. Sample silicon concentration fields for growth of dilute SiGe alloy (Rag = 4 pm/s.

=

0) in Grenoble furnace with the graphite ampoule;

P.M. Adornato, R.A. Brown

RaT’ I x Ra S ‘—Ix

106 I0~

/

RaT’ I xI06 Ros~’IxI06

I H

183

Convection and segregation in directional solidification

Pill x

ROT’ I

06

Ra 5’-I x I0~

_

6

.~rx

-9.78x10

I0~~Io-5 &

Fig. 25. Flow fields for growth of SiGe alloy in Grenoble furnace with graphite ampoule; RaT = i x i06.

ROT=1x106

Ras ‘—lx tO

~

I 0.662

Ra-Ix’~6 T

0.666

R0

5”I

X

106 I

_____

_____

______

______

H

H

6 ROT-1x10 ROg’H X I0~

0.5701

0.4246

O.2792~ _

4-

II 4.

4..

Fig. 26. Silicon concentration fields for growth of SiGe alloy in Grenoble furnace with graphite ampoule; RaT = 1 x 106.

184

P.M. A dornato, R.A. Brown

I

1=

Ra1’ lx 10

~max

x

I

(~

Ix I0~

Ras ~Ix

106

ROTrIXIO7 Ra5’Ix I0~

(H

(1I

I

Ras’lx I0~

‘.35

Convection and segregation in directional solidification

RO

~, I

/

=

~

I

~

~max’

1.89 x

~

~i’mln= -9,BQx 0 ~==fl~~j

‘k(T~Ifl

-6.82 x I0—1~

max

1

‘Pmin =

-1.90 x I0~

1

t

4-

Fig. 27. Flow fields for growth of SiGe alloy in Grenoble furnace with graphite ampoule; RaT = 1

H

_

RaT’ I x I0~ R0sIX

I0~

ROT:I L0.52I0~

r ~

0.3915

I

x I0~

Ras

-

x iO.

_

_____

6

-—

Ix 10

0.54c~

Ra T

I x I0~

0.5695

Ras ~Ix I0~

______

I

0.3995 0.42 38 0.26 0.2585 0.2780

Ii

—

1.88 x I0~

‘—

4-

Fig. 28. Silicon concentration fields for growth of SiGe alloy in Grenoble furnace with graphite ampoule; RaT = 1

x

‘°

P.M. Aaornato, R.A. Brown

/ Convection and segregation

meshes, we cannot discriminate between a physical loss of axisymmetric, steady-state solution and error in the finite element approximation as the cause of this failure. The silicon concentration fields which correspond to the flows shown in fig. 25 are given in fig. 26. The small degree of convective mixing apparent for the lowest value of Ra~is lost as the silicon concentration is increased. The onset of the axially-stacked cells at the highest value of Ra8 is not represented in the concentration fields near the interface, and thus can not be detected by examining either the radial or axial solute distributions in the crystal. In this sense, the appearance of the transition to the axially-stacked cells has no effect on the crystal growth. Adding silicon to a system at a higher thermal Rayleigh number shows the same transition to the axial cells caused by the instability of Hart; streamline for RaT = 1 iO~and three values of Ra are shown in fig. 27. The concentration fields corresponding to these flows are shown as fig. 28. The core of well-mixed melt adjacent to the melt/crystal interface for Ra5 = 0 shrinks with decreasing Ra until at Ra = 1 x iO~the concentration field resembles the one appropriate for diffusion-controlled growth. Again, the onset of the additional cells does not change the concentraX

—

—~

j

0.2C I

I

I

—

V9 ‘4pm /8

..~.

0.25

.

~

0.24

‘

0.23

.

0

0.22

m

.

~—m~O.I6I 0.21

AMPOULE WALL

0

0.02

0.04

0.06

~

DIMENSIONLESS RADIAL COORDINATE r

Fig. 29. Effect of hquidus slope on rnterfacial solute concentration profile c( r, h ( r)) for GeSi alloy growth in Grenoble furnace; curves are for m = 0 and m = 0.161, the value appropriate for 5% Si in Ge.

in directional solidification

185

tion field near the interface. As shown in fig. 29, increasing the dimensionless slope of the liquidus curve m from zero to 0.161, the value appropriate for 5% Si of Ge, leads to an increase in the deflection of the interfacial concentration c(r, h(r)), but had little other effect on the results. Further increases in m eventually lead to a numerical instability that cause spatial oscillations in the interface shape at the same frequency as the mesh. The relationship between this instability and constitutional supercooling is discussed in ref. [20].

7. Comparison with experimental measurements in MIT system A detailed comparison between the calculations reported here and the experimental measurements of Wang and Witt [14] was carried out to quantify the validity of our model. The calculations were compared directly to results of crystal growth experiments for (a) the shape of the melt/crystal interface measured by Peltier demarcation, (b) the radial segregation of gallium across the ampoule, and (c) the degree of axial segregation measured by the effective segregation coefficient as a function of distance grown. The measurements of each of these properties are given in ref. [35] and only the comparison with our calculation is reported here. The thermal boundary conditions used to model the ampoule, furnace and their interactions were confirmed by comparing calculations with temperature measurements for a solid graphite rod placed in the furnace to simulate an ampoule and crystal. The predicted and measured temperatures are shown asfunction fig. 30 for the furnace profile and number described in section 2.3. Biot The average error in the predicted temperatures was only 0.011 in dimensionless units, without any attempt at fitting the heat transfer coefficients to improve the agreement. Direct measurements of temperatures in the actual GaGe/ampoule system is difficult, so companson of the thermal field there to the computed value relied on the prediction of the melt/crystal interface location. Interface shapes predicted by .

186

P.M. Adornato, R.A. Brown 3 0

I

I

I

/

Convection and segregation in directional solidification

I

—

_.

—

•

EXPERIMENTAL MEASUREMENTS

—

FINITE ELEMENT CALCULATIONS

1.5—

~

~

-

EXPERIMENTAL

~FINITE

—

I

U

ELEMENT

I380~m

-

0.00

CALCULATION

I

I

I

I

0.02

0.04

0.06

0.08

0.10

DIMENSIONLESS RADIUS, 0

-I

~

—

-

—

—

0.2

0.4

0.6

o~o

1.0

DIMENSIONLESS TEMPERATURE 8(o,z) Fig. 30. Comparison of axial temperature profiles along the center of a graphite rod predicted by fimte element analysis and measured by Wang [35]

the calculations and measured in the experiment are compared in fig. 31. Again, the agreement is good, considering that heat transfer parameters have not been adjusted to fit the data, Prediction of the axial and radial concentration profiles in the crystal growth experiment are the most critical comparison which can be made between experiment and calculation. The Grashof number, defined by eq. (19b), was calculated as Gr 1.86 x iO~ for the experimental system of Wang [35]. The ampoule used in the measurements had length 30 cm with the other dimensions of the system given in table 3. The calculations presented in sections 5.1 and 5.2. indicate that most of the additional length of melt beyond the shorter value used in those calculations has not influence on the solute field, because the gallium concentration is essentially constant above the junction of the adiabatic and hot zones. We expect

Fig. 31. Comparison of melt/solid interface shapes in growth of GaGe in MIT furnace predicted by finite element analysis and measured by Wang [35] for two sets of hot and cold zone temperatures.

that the upper flow cell will fill a longer ampoule and continue to mix the solute at this high value of Gr. We compare our calculations with the shorter ampoule directly to Wang’s experiments. Unfortunately, the Petrov—Galerkin formulation did not yield accurate results for the concentration field at the thermal Rayleigh number RaT = 3 x 10~,which is appropriate for modeling the experiments of Wang and Witt, with L = 7.62 cm. The differences between the magnitude of the radial segregation ~c predicted by the calculations and the measurements of Wang and Witt [14] are shown in fig. 14 for the three crystal growth rates of 2, 4, and 8 fim/s. Notice that the differences between the amount of segregation predicted for each growth rate by the calculation are in the same proportion as in the experiments. Also, the comparison between the experiments and calculations shows that the segregation in the experiments is indicative of intense convective mixing, so that differences in concentration are confined to thin boundary layers along the melt/crystal interface and to internal layers separating the axially stacked convection cells. We have used the apparent boundary layer structure of the concentration field to predict the

P.M. A dornato, R.A. Brown

/ Convection and segregation

radial segregation for higher values of RaT. To do this we have assumed that the concentration field along the interface within the boundary-layer will vary according to the scaling law

187

in directional solidification

modelled by the expression for RaT>> 1 &

(75.5J’

where

+

10)(RaT)~”7,

(33)

has units of ~.Lm/s.

c 1(r) — —

C(r) (Pe)’~ C(r) (Sc/Pr)” (Gr)” ~ ‘R ~ 0~r, aT) ~~32) ‘~

,.

~

‘

where the dimensionless groups are defined in table 2. The radial concentration difference & also is assumed to scale as ~1c= (&)~~ (Ra)” Table 5 lists the values of n and (&)~computed from the finite-element calculations for the largest two values of RaT and the resulting value of the

Axial segregation is predicted from the calculations by computing the effective segregation coefficient according to the formula eq. (1). Three results also are extrapolated to the convection level for the MIT furnace by using the correlation (34)

keff = Keff(RaT)nk.

The extrapolated values of keff are given in table 5. The values of k~ffreported by Wang [35]were

radial segregation obtained by extrapolating this

measured from a single crystal growth experiment,

result to RaT = 3 x 10~for the three growth rates. These values are higher than the experimental results for all three growth rates, but give confidence that the structure of the flow and solute fields adjacent to the interface are accurately predicted by the finite element calculations. At these high thermal Rayleigh numbers the radial segregation is most appreciable near the ampoule wall. If spreading resistance measurements were not available near the outer radius of the crystal, the reported experimental values would be lower. The values of n computed for the three values of the growth rate varied only between (—0.137 and 0.148 indicating that the structure of segregation is relatively insensitive to changes in J’. The pre-exponential constant for eq. (32) varied linearly with V~,as expected from fig. 14. Then the’ radial segregation in our GaGe calculations can be

in which the length of the melt decreased by a factor of two between the measurements. The range in values of kett shows the dependence of the flow and solute fields on the length of the melt. Wang found that the measured values of keff were in good agreement with the normal freezing equation [21] which assumes complete mixing in the melt. Changing the length of the ampoule in the finite element calculations would affect the value of keff. The added length would increase the size of the upper flow cell which has the lower core concentration (see fig. 9), thereby decreasing ((c)) and increasing keff. The effective segregation coefficient should be independent of melt length when the mixing is intense, so this trend is not correct. Axial segregation in the MIT furnace is controlled by the composition in the core of the lower

—

Table 5 Coefficients obtained 8 experimental by fitting values GaGe of ~1c calculations and k~ in MIT furnace to eqs. (32) and (35) and extrapolating values of & and k~~fto RaT = 3 x i0 1~ 7 from ref. 1351 are also listed; the range of thermal Rayleigh numbers used in the fitting of the equations was ix ~ ~ Ra.~.~ 5 X i0 Growth/rate J’~ Parameter Constant Exponent Extrapolation Experimental (pm/s) InC 8 value[35] 0 to3xlO (%) (%) 2

.~c

~ 4

i.Sc

8

z~c

k~’~f

5.08 —1.28 5.62 —0.404 6.42

—0.148 —0.057 —0.137 —0.085 —0.141

8.95 0.092 19.18 0.128 39.06

0.697

—0.123

0.181

4 0.090—0.095 14 0.104—0.111 30

k~ 1

—

P.M.

188

Adornato, R.A. Brown

/ Convection

cell near the interface, not the composition in the hot zone. This suggests the alternative definition of the effective segregation coefficient for the PSSM as elf

k(c) tin /c*

‘

(35~ ‘

where c” is the concentration in the cell adjacent to the interface. Extrapolating the value of k~,to RaT = 3 X iO~according to ta scaling law similar to eq. (35) gives the coefficients and values listed in table 5. These values agree well with the experiment, reinforcing the idea that mixing in the heat-pipe furnace is limited by transport across the internal layer between the cells, and that the composition is not uniform throughout the melt, even with intense laminar convection.

8. Discussion Thermosolutal convection in directional solidification arises from the interaction of driving forces determined by both the alloy’s thermophysical properties and by the furnace and ampoule design. The flow pattern and concentration fields have extremely rich structures. More importantly, the details of this structure sets the axial and radial segregation of impurities in the crystal. The vertical Bridgman system used at MIT leads to a two-cell flow structure; with one cell driven by the radial gradients at the junction of the adiabatic and hot zones of the furnace and another adjacent to the interface caused by the differences in thermal conductivities of melt, crystal, and ampoule. The constant gradient furnace designed by the Grenoble group has only a single cell near the interface. In the solidification system studied here, solute is transported mostly by diffusion between adjacent cells, because the convective motion is laminar, Then radial segregation is set primarily by the flow adjacent to the interface when this cell is long enough to completely contain the solute gradient for diffusion-controlled growth. The upper cell in the MIT furnace plays little role in determining the solute distribution along the interface. Axial segregation depends on the intensity of convective mixing throughout the ampoule, so that the cellu-

and segregation in directional solidification

lar structure of the flow pattern is important. Solute transport between flow cells in the heat-pipe system limits the effectiveness of this mixing and leads to high solute concentrations in the cell next to interface. Thecell, degree axial segregation is set the mainly by this as of demonstrated by the comparison between the effective segregation coefficients calculated here and measured by Wang [35]. The thermal conductivity of the ampoule material and the thickness of the ampoule wall can be used as parameters to change the structure of the flow, as well as to modulate the flow intensity. For example, the shape of the single buoyancy-driven cell in the constant gradient furnace is influenced by the choice of ampoule material. For the graphite ampoule, the radial temperature gradients near the interface are large and convection is mostly confined to the melt adjacent to the interface. For the composite ampoule, the radial temperature gradients near the interface are smaller and convection is more uniform along the ampoule wall and the flow cell extends further from the melt/crystal interface along the ampoule wall. We demonstrated the changes in the flow intensity in the constant gradient furnace by calculations with the graphite and composite ampoules and in the heat-pipe geometry by the calculations with the boron nitride and quartz ampoules. The dependence of the degree of radial segregation on the thermal Rayleigh number RaT in each of the system studied here can be characterized by the curves shown in fig. 20. Segregation is lowest for either diffusion-controlled growth (small RaT) or when the intense laminar mixing is present (large RaT), but reaches a maximum for intermediate values of flow intensity. Changes in the ampoule material, in the furnace design, and in the addition of a stabilizing solute all alter the values of segregation along this curve, but do not change its qualitative form or interpretation. This result is crucial for understanding the effects on solute segregation.of microgravity sohidification and of applied magnetic fields. In both situations the body forces acting on the melt are altered in an attempt to damp convection. Both approaches are clearly effective at eliminating

P.M. Adornato, R.A. Brown

/

Convection and segregation in directional solidification

time-periodic and chaotic convection present in larger-sized systems, but the effect of the reduced convection level on solute segregation has not been well understood. Our calculations clarify this point. If the reduced gravitational field has constant magnitude and is aligned with the crystal axis, it is modelled exactly by reducing the thermal Rayleigh number. The appropriate value for the experimental system of Wang [35] is Ra~= (g*/g)Ra~, where the asterisk denotes variables appropriate for the reduced gravity environment. It is clear from examining fig. 20 that a ratio of (g*/g) of 10_2~_10_3 may lead to more radial segregation in the crystal than was present in the earth-bound experiment. The decrease is gravitational level needed to reach diffusion-controlled growth increases with R3, so that the effect is more severe for larger ampoules. The effect on radial segregation of the damping of convection by an axially-aligned magnetic field can be visualized in terms of a plot analogous to fig. 20 but with the abscissa RaT replaced by the reciprocal of the dimensionless strength of the field. Curves of constant RaT, e.g. crystal growth on earth, will have the maxima characteristic of the curves in fig. 20, so that intermediate values of the field strength will lead to poor radial uniformity. Calculations demonstrating this feature for a realistic crystal growth system will be presented later.

Acknowledgements

This research was supported by the Microgravity Sciences and Space Processing Program of the United States National Aeronautics and Space Administration, and by a grant from the Department of Energy for use of the computer facilities at Los Alamos National Scientific Laboratory. We are also indebted to C.A. Wang and J.J. Favier for many valuable discussions concerning directional solidification of semiconductor materials and the designs of their crystal growth systems. The cornments of D.H. Kim on an earlier version of this paper also are appreciated.

189

References [1] CE. Chang and W.R. Wilcox, J. Crystal Growth 21(1974) 135. [2] P.C. Sukanek, J. Crystal Growth 58 (1982) 208. [3] P.C. Sukanek, J. Crystal Growth 58 (1982) 219. [4] T.J. Jasinski, W.M. Rohsenow and A.F. Witt, J Crystal Growth 61(1983) 339. [5] T.-W. Fu and W.R. Wilcox, J. Crystal Growth 48 (1980) 416. [6] R.J. Naumann, J. Crystal Growth 58 (1982) 589. [7] R.J. Naumann and S.L. Lehoczky, J. Crystal Growth 61 (1983) 707 [8] T.J. Jasinski, W.M. Rohsenow and A.F. Witt, J. Crystal Growth 71 (1985) 295. [9] C.J. Chang and R.A. Brown, J. Crystal Growth 63 (1983) 343. [10] C.J. Chang and R.A. Brown, J. Computational Phys. 53 (1984)Brown, 1. [11] R.A. C.J. Chang and P.M. Adornato, in: Model!ing of Casting and Welding Processes, Eds. J.A. Dantzig and J.T. Berry, AIME, 1984. [12] F.M. Carlson, AL. Fripp and R.K. Crouch, J. Crystal Growth 68 (1984) 747. [13] G B. McFadden, R.G. Rehm, SR. Coriell, W. Clark and K.A. Morrish, Metall. Trans. A15 (1984) 2125. [14] CA. Wang and A.F. Witt, Annual Report Materials Processing Center, Massachusetts Institute of Technology, 1984. [15] A. Rouzaud, D. Camel and J.J. Favier, J. Crystal Growth, to be published. [16] CA. Wang and A.F. Witt, JR. Carruthers, J. Crystal Growth 66 (1984) 299. [17] ED. Bourret, J.J. Derby and R.A. Brown. J. Crystal Growth 71 (1985) 587. [18] I. Christie, D.F. Griffiths, AR. Mitchell and O.C. Zienkiewicz, Intern. J. Numer. Methods Eng. 10 (1976) 1389. [19] AN. Brooks and T.J.R. Hughes, Computer Methods AppI. Mech. Eng. 32 (1982) 199. [20] P.M. Adornato and R.A. Brown, Intern. J. Numer. Methods Fluids, in press. [21] E. Scheil, Z. Metallk. 34 (1942) 70. [22] PB. Rhines and W R Young, J. Fluid Mech. 133 (1983) 133. [23] C. Huh and L.E. Scriven, J. Colloid Interface Sci. 35 (1971) 85. [24] T. E.B.Jasinski, Dussan,PhD Ann.Thesis, Rev. Fluid Mech. 11(1979)371. [25] Massachusetts Institute of Technology (1982). [26] SM. Pimputkar and S. Ostrach, J. Crystal Growth 55 (1981) 614. [27] D. Camel and J.J. Favier, J. Crystal Growth 67 (1984) 42. [28] A. Rouzaud, D. Camel and J.J. Favier, J. Crystal Growth 73 (1985) 149. [29] G.M. Harriott and R.A. Brown, J. Fluid Mech. 126 (1983) 269. [30]R.F. Dressler, J. Crystal Growth 54 (1981) 523.

190

P.M. Adornato, R.A. Brown

/ Convection and segregation

[31] Y.F. Pan and A. Acrivos, Intern. J. Heat Mass Transfer 11 (1968) 439. [32] A. Acrivos, Chem. Eng. Sci. 21(1966)343. [33] H. Schlichting, Boundary Layer Theory (McGraw-Hill, New York, 1955). [34]S. Kimura and A. Bejan, Phys. Fluids 28 (1985) 2980. [35] C.A. Wang, PhD Thesis, Massachusetts Institute of Technology (1984). [36]P. Huyakorn, C. Taylor, R. Lee and P. Gresho, Cornputers Fluids 6 (1978) 387. [37]H.M. Ettouney and R.A. Brown, J. Computational Phys. 49 (1983) 118. [38] O.C. Zienkiewicz and J.C. Heinrich, in: Finite Elements in Fluids, Vol. 3, Eds. Gallagher, Zienkiewicz, Oden, Morandi-Cecchi and Taylor (Wiley, New York, 1978).

in directional solidification

[39] J.J. Derby and R.A. Brown, Chem. Eng Sci. 41(1986) 37. [40] L.H. Ungar and R.A. Brown, J. Computational Phys., submitted. [41] P. Hood, mt. Intern. J. Numerical Methods Eng. 10 (1976) 379. [42]J.E. Hart, J. Fluid Mech. 49 (1971) 279. [43] RB. Bird, WE. Stewart and EN. Lightfoot, Transport Phenomena (Wiley, New York, 1960). [44]CF. Chen and D.H. Johnson, J. Fluid Mech. 138 (1984) 405. [45]R.K. Crouch, AL. Fripp, W.J. Debnam, RE. Taylor and H. Groot, presented at Materials Research Society Meeting, Boston, November 1981.