- Email: [email protected]
Diffusive-convective vapor transport across horizontal and inclined rectangular enclosures
Journal of Crystal Growth 67 (1984) 241—254 North-Holland, Amsterdam
241
DIFFUSIVE-CONVECTIVE VAPOR TRANSPORT ACROSS HORIZONTAL AND INCLINED RECTANGULAR ENCLOSURES B.L. MARKHAM and F. ROSENBERGER Department of Physics, University of Utah, Salt Lake City, Utah 84112, USA
Received 23 August 1983; manuscript received in final form 13 February 1984
The enhancement of vapor transport across horizontal and inclined rectangular (2D) enclosures by expansive convection is studied numerically for a range in Grashof number that extends up to boundary layer flow. Comparison is made with predictions of the simplifying model and resulting analytical treatment of Klosse and Ullersma for various aspect ratios of the enclosure. The influence of varying Schmidt and Prandtl numbers is investigated. Solutions for interfacial (growth) flux and temperature distributions are given, showing practically important non-uniformities.
I. Introduction There has been considerable interest in accurate predictions of transport rates for crystal growth from vapors in closed ampoules. In particular, this interest has grown since chemical vapor transport experiments recently conducted under low-gravity conditions aboard space-crafts gave results [1,2] that still defy quantitative interpretation. The earlier, mostly 1D or idealized 2D treatments of vapor transport (VT) in closed ampoules have been reviewed elsewhere [3,4]. Among these, the most realistic and, hence, widely used model for diffusive—convective horizontal VT is due to Kiosse and Ullersma [5] (hereafter referred to as KU). Various experimental studies, however, failed to quantitatively confirm KU’s predictions for a significant range of characteristic parameters; see, e.g., refs. [6,7]. Most of the simplifying assumptions made by KU to obtain a closed form analyti cal solution were relaxed in later numerical work [8]. Due to numerical difficulties, however, the Boussinesq approximation was used and only flows with rather low Grashof number were considered. For this limited parameter range good agreement with the simpler KU model was obtained [8]. In the present paper, based on an improved numerical approach, we cover the whole convection
(Grashof number) range that is of interest in the laboratory. In addition to the horizontal transport problem we consider the inclination of the rectangular enclosure with respect to gravity over the whole range between “top-heating” and “bottomheating”. Stimulated by some of the simplifications in KU’s work, we have also explored the influence on transport rates of variations in vapor properties (Schmidt and Prandtl numbers) over a range that is representative of vapor mixtures in the VT practice. Since KU’s work forms the starting point for our efforts we shall give a brief summary of their model and approach for solving it. A rectangular horizontal cavity of height h and horizontal extent L is assumed, see fig. 1. The third dimension is taken as very large. The vertical end walls are ___________
L
-___________________________
h T5
T C
Fig. 1. Expansive convection due to horizontal temperature gradient in rectangular enclosure. Overall flow, and velocity profile (S-profile) in vertical Center plane; schematic.
0022-0248/84/$03.OO © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
242
B, 1.. Markham. f Rosenherger / l)s,/,fuses ‘e --~‘ons’ee1,s’elapor transport
differentially heated to 1~and J~.and drive expansive convection of the vapors between them. The end walls form also the interfaces between solid source and vapor (at L,), and vapor and crystal (at 7~.),respectively. In solving for the mass transport rate, KU disregard the heat transfer across the enclosure. A linear temperature gradient aT/ax = (7’~— 7’~)/Lis assumed in the vapor, irrespective of the actual convection conditions. Note that with this assumption it is implied that the vapor has infinite thermal diffusivity and, hence, vanishing Prandtl number, Pr (see nomenclature at end of paper). Solutal density gradients are ignored. Hence, the fixed temperature gradient results also in a fixed density gradient ap/ax. Furthermore, velocity contributions arising from the net mass transfer (vapor in-flux and out-flux at source and crystal, respectively) are also ignored. The convective flow is then obtained by KU by assuming a stream function for an inifitely long, horizontal channel with anti-parallel flow (see the S-profile in fig. 1). The boundary conditions for the vertical end walls were grafted onto the stream function and the resulting unknown parameters were solved for with a variational technique. This newly found velocity distribution was then used to solve for the concentration profile. The integral solution allowed KU to make predictions about the enhancement of mass transport due to convection in the form of a rather simple mathematical function of aspect ratio AR = L/h. Grashof number Gr and Schmidt number Sch. No solutions for local concentration distribution and growth flux were given, Intuitively one can expect the following limitations from these model assumptions and the method of solution. Since a specific symmetric stream function was chosen, the model cannot account for deformations of the velocity field (development of boundary layer behavior) that arise in reality at higher flow velocities i e higher Or Synonymous with this restriction is the assumption of a constant density gradient (equivalent to Pr = 0). This restricts the drive for convection to the core mechanism. In reality, T(x) deviates increasingly from the assumed linear profile with increasing convection strength, and eventually boundary layer entrainment provides the dominant drive at
high Or (9J. For vapors, where Pr is of order one. the restriction to the core mechanism can he cxpected to he a serious limitation even at not too strong convection conditions. Also, since the stream function assumed for the core has the character of the flow in an infinitely long tube, one can expect poor results even at modest Or for systems of small aspect ratio. In the following we will quantitatively explore these questions. together with an expansion of the model to arbitrary orientations with respect to gravity, and present solutions for concentration, temperature and growth flux distributions.
2. Model and theoretical analysis We shall consider a rectangular cavity of height h and length L, the end walls of which are inclined with angle ~ + s~/2 with respect to the gravity vector, see fig. 2. The cavity is assumed sufficiently deep so that the third dimension may he neglected. As in the KU model (fig. 1) transport takes place between two flat, isothermal interfaces that are kept at 7’~,and 7., respectively, with 1~>i~.(‘omponent A sublirnes from the source, is transported through and with a non-reactive component B, and condenses at the crystal interface. We assume complete rejection of B and equilibrium of A at both interfaces. The bounding walls are assumed to be impermeable to both components and thermally ideally conducting. Radiative heat transfer is neglected. The molecular weights of the two vapor components are assumed equal, MA = M 15, which restricts our model to expansive (“thermal”) convection. Also, frictional contributions to energy
~ Source~
1
I
j
——_____
Vapors
_______
~~taI
~
~B [x
/ Fig. 2 Model for diffusive—convective vapor transport across rectangular enclosure of arbitrary orientation ~ with respect to gravity.
B. L. Markham, F. Rosenberger / Diffusive — convective vapor transport
transport are ignored. The viscosity, binary diffusivity and thermal diffusivity are assumed concentration and temperature independent, The transport in this cavity is governed by the system of nonlinear, coupled conservation equations for momentum (Navier—Stokes), mass (continuity), species (diffusion) and energy. These equations are, respectively, in rectangular coordinates and steady state
a
a
~(puu)
=
+~(pvu)
a~
F
a
au
2 /
au
av il
a2v + ~ a2u ayax —~+pg~,
ay
ay
—
rigidity v(x,0)
=
v(x,h)
=
0,
(6)
no-slip u(x,0)
=
u(x,h)
=
0,
(7)
conducting
(Ia)
a~
av
2 /
au
av
a
(lb)
0,
(2)
—
T,,—T~ L
(8)
a
au’ç
1
—
—
DAB awA(o,Y) W,~(0,y) ax
‘
(9)
v(0,y)=0,
(10)
uniform B-concentration (11)
uniform temperature
a
awç \
/
~
ay a(pvT)
=
pB(0,y)=const.,
-~—(puW,)+-~—(pvWç)
~(puT)+ a
—
—
no-slip
a2v ax
=
u(0,y)
~
+j.t—+j.t——~+pg5’,
a(pu) a(pv) ax + ay
aT(x,h) ax
—
out-diffusion of A and impermeability to B
a F
a2u axay
aT(x,o) ax Source:
a(puv) a(puv) ax + ay =
the transport region. The boundary conditions are: Walls: impermeability awA(X,o) — 8J’VA(x,h) = 0, (5)
ax +1L~
243
2T
~
k (a
a2T
T(0,y)
=
T 5.
(12)
(3)
Crystal:
(4)
in-diffusion of A and to B DABimpermeability aw~(L,y) =
-
u(L,y) As in our work on diffusive—convective transport in vertical, cylindrical ampoules [3,10], we have
kept the “compressible terms” in the momentum transport equations. Also, neither creeping flow nor boundary layer assumptions were made, in order to retain applicability to the wide range in flow parameters (Gr) to be investigated. In addition to (1)—(4), the ideal gas law and Dalton’s law of partial pressures were used to link density to temperature and concentration. To solve these equations we must state boundary conditions on both walls and interfaces enclosing
no-slip v(L,y)
1— WA(L,y)
=
0,
ax
(13)
(14)
uniform B-concentration pB(L,y) = const.,
(is)
uniform temperature T( L y) T~.
(16)
,
=
Clearly, some of these boundary conditions merit more discussion. Perfectly conducting walls
244
B.L. Markham, F Rasenherger
(or, for that matter, perfectly insulating walls) cannot be obtained in a crystal growth process. Hence, (8) must be understood as one approximation which is probably reasonable, comparing, for instance, the thermal conductivities of quartz and typical gases. The no-slip condition (7), has been
/
Diffusive
-
conteclIse taper transport
molecular weight of the vapor components any diffusion slip can safely be excluded [111. With MA = MB we can also neglect solutal convection effects. As we have pointed out before [3,12] there is a singularity at each corner of the cavity. The singularity arises from the simultaneous statement of the no-slip condition and of in- and out-diffusion, respectively, at these points. For a thorough discussion of this problem the reader is directed to ref. [12]. Suffice it to say that these singularities
out to he highly three-dimensional, even in horizontal cylinders with ideal linear temperature profile. ‘This forms probably the strongest limitation of the following results and discussion. However, the tools (codes and computer power) for a full-fledged 3D treatment of this problem including momentum, energy and species transport. are Just becoming available. Hence, the semi-quantitative guidance that can he obtained from our resuits, can probably he useful for the practitioner until more rigorous treatments can be carried out for the high Or conditions often found in the laboratory. The system of non-linear, coupled partial differential equations (1)--(4) was solved with a finite difference method. We developed an algorithm which is a variation of the code TEACH [171. The grid system used had 16 and 32 grid lines in the r-
appear to have little effect on the solutions given here. In reality the interface shape is a result of the A-flux distribution, i.e. for non-uniform interfacial
and x-direction, respectively. For some of the soiutions, in particular at higher Or, we used non-uniform grid spacings in the x-direction. with tighter grid packing close to the interfaces, where higher
flux, the shape of the interface will change and not remain flat as we have assumed. The evolution of the interface shape as a function of the prevailing
gradients in velocity, concentration and temperature are encountered. As starting solutions for the iteration process we used typically solutions to the corresponding I D problem with linear temperature gradient. However, for high Or conditions, this approach offers no advantage over an “initial fill” with zeroes and a “boot-strapping” procedure from lower to higher Gr becomes very important for efficient computation. During a computational run the total pressure was held constant only at one point, namely at the intersection of the third x- and y-grid lines (safely spaced from the corner singularity). The pressure at all other grid points was up-dated with each iteration to satisfy continuity. The criteria for convergence were an unchanging velocity field and balance of the mass flux for all cross-sections x = const. of the cavity. The calculations were carried out on a VAX 750 computer. Double precision was used to overcome certain truncation instabilities. The results presented here, including extensive program development runs, required approximately 300 CPU hours. The establishment of the physical integrity of any numerical scheme is a non-trivial task. This is particularly true for the systems considered here due to a lack of non-intrusive measurements of
the subject of some discussion, as summarized in ref. [3]. However, with the assumption of equal
heat and mass transfer conditions (which, in turn, change with the interface shape) forms a formidable mathematical task (Stefan problem). However, the rate at which the interface moves is typically more than three orders of magnitude slower than representative vapor flow velocities. Hence, the assumption of an unchanging shape of source and crystal is reasonable. To assume planar interfaces at the respective ends of the ampoule, however, is not necessarily representative of the crystal growth practice. Most workers obtain many little crystals in these ampoule regions, unless specific provisions (seeding, compacted source, etc.) for planar interfaces are made. Yet, since the purpose of this paper consists in pointing out useful trends, we feel comfortable with this mathematically greatly simplifying assumption. During the course of this work, due to input from our parallel experimental efforts [13,14] and from 3D calculations [15,16] it has become obvious that, in general, 2D treatments of VT processes in cylindrical ampoules must be rather limited. In particular at higher Gr, the convective flow turns
B. L. Markham, F. Rosenberger
/
temperature and concentration distributions in closed vapor transport ampoules. Velocity distributions have become available only recently
40
of the experimentally obtained convection behavior. These efforts, together with the results of other numerical schemes, will be reported elsewhere [16]. For the concentration distribution obtained with our scheme we gained confidence from other runs that we tailored to simulate earlier diffusive—convective experiments with gases in semi-open systems [18]. This shall be discussed in connection with fig. 17. In the following we will present the results grouped according to effects of (a) changing temperature differences and pressure, i.e. changing Grashof number; (b) changing materials properties of the vapor, i.e. changing Prandtl and Schmidt numbers; (c) changing geometry, i.e. changing ampoule length, and orientation with respect to gravity.
‘‘‘i
‘‘‘‘‘~‘‘“
Pr
[13,14]. Hence, ir order to test our scheme, we have conducted several numerical runs with single component fills, i.e. with component B only. Good
agreement was obtained with certain 2D features
245
Diffusive — convective vapor transport
30
8
= 5
Pe
~
i.56
20
0 7
1
Sch AR
-
~
—
-
•
-
-
~ ~ -
-
KU 1
-
‘‘‘‘‘13
‘‘‘‘~~~
10
I ~‘~‘‘‘‘i’~~
I
-~
106
GRASHOF NUMBER
Fig. 3. Variation of Sherwood number with Grashof number for several Pe’s at AR = 5. Heavy line from our model calculation, lighter line prediction from KU model.
Fig. 3 illustrates, for AR = 5, the importance of convection for mass transport in VT. For Or ~ 2 x i0~, diffusion dominates, as reflected by
3. Results and discussion
60
The following dimensionless groups are stated with each result to define the system: Grashof number Gr, Prandtl number Pr, Schmidt number
50
Sch, Péclet number Pc, aspect ratio AR. In addition the angle 4) (see fig. 2) is given. The density, that enters Or, Pr and Sch, was taken at the pressure monitoring point for the average temperature(1,+7~)/2.
40
-
-
30
-
-
3.1. Effects of Grashof number
20
-
‘I
-
‘‘‘‘I
‘‘‘‘‘I’’’!’
~
-
AR
=
KU
These solutions are all for “horizontal ampoules”, i.e. for 4) 90°.Aspect ratios of 5 and =
10 are considered. Constant values for Pr
=
and Sch 1.0 were used throughout this section. Most of the Gr range was scanned by changing the (monitor point) pressure. A few low Gr results were obtained by increasing the viscosity. For comparison, a representative experimental range in Gr is from below i04—2 X 10~’[5,6]. =
10
-
-
0.7 0
.1
io~
,,,,
I
106
GRASHOF NUMBER Fig 4. Variation of Sherwood number with Grashof number for Pe 0.95 at AR 10. Heavy line: our model, lighter line: KU.
246
B. L. Markham, F Rosenherger / Diffusive
Sherwood numbers, Sh close to unity. For Or> 2 >< iO~, however, the convective contribution be-
comes dominant. For more slender “ampoules”, e.g. for AR = 10, as shown in fig. 4, the cross-over from the diffusion- to the convection-dominated regime occurs at a somewhat higher Gr. This is intuititively obvious, since with increasing AR there is viscous interaction of the flow with a longer sidewall, The corresponding KU solutions are also plotted in figs. 3 and 4. One sees that at Gr’s that are in the middle of the experimental range, KU’s predictions level out, whereas our less restricted model yields a continuous increase of Sh with Or, in good agreement with experimental observations [7]. The steep section of our curves (Or ~ 5 >< 10~) represent functions Sh ~ Gr°, where the power n is 0.28 and 0.25 for figs. 3 and 4, respectively. This compares very well with heat transfer results obtamed from boundary layer behavior [9,19]. We have also investigated the influence of the Péclet number on the Sh, i.e. the influence of the in-flow, and out-flow velocities from net mass transfer. A derivation of the relation for the Pc used here (see Nomenclature) is given in the appendix of ref. [3]. Somewhat surprising to us, a variation of 0.95 ~ Pc ~ 0.56 (which well covers the experimental range) has practically no result on Sh (fig. 3). This reflects the negligible effect of the diffusive through-flow on total velocity even in the presence of only weak convection, The corresponding solutions for the velocity field show that even at Gr = I0~,where convection aids net transport very little (i.e. Sh 1), the convective velocities are two orders of magnitude larger than the diffusive velocities. Because of this insensitivity of the results to reasonable Pc variations we have used a constant Pc = 0.95 for all following caiculaAt this point one may be tempted to neglect the species diffusion equation altogether during solving for velocity and temperature distributions and, then to determine the concentration distribution based on these uncoupled solutions. This would reduce the computation time significantly. One could even get fairly reliable information on the local fluxes at the interfaces in this way. However,
this procedure neglects whatever effect these in-
convective vapor transport
—
and out-fluxes might have on the overall flow. As we shall see, the interfacial fluxes can change the entire form of the flow. Hence, the species diffusion equation must remain coupled during the solution process. Let us now concentrate on the physical reasons for the unrealistic results of the KU approach. Fig. 5 shows computed temperature profiles taken at half height (v = h/2) along the length of the cayity, for AR = 5 and three different Or’s. At Or = 5 x 1O~,where we find excellent agreement with the Sh predicted by KU, the temperature profile is practically linear. At Or = 1.78 x iO~,where the transport results of the two models begin to significantly diverge, the temperature profile begins to he deformed by convection. At Or = 1,14 X l0~, representing the highest valtie in fig. 3, and the largest disparity with KU’s prediction. the temperature profile has assumed boundary-layer type shape. Such strong deviations from the linear distribution must, of course, lead to failure of the KU model. Whereas in reality (and in our model) a transition from core-driven to boundary-layerdriven flow occurs, the KU model is by assumption limited to core-driven flow at all Or’s. Additional insight can be obtained from computed velocity profiles. Fig. 6 represents core flow velocities u(t’) at x = L/2 for the Or values of fig. 5. We see that only for the lowest Or one obtains a symmetric “S-profile”. At higher Or the cross-over from positive to negative velocity occurs signifi-
T
F
‘~0.8
~ ~
o.~
Or
~
.14~ -~
-
~ Pr ~
0.2
0
0
0 7
SeS
=-
1.0
AR
=
5 0.2
I
I 0.4
0.6
0.8
1.0
x/L
Fig. 5. Temperature profiles at h /2 for Gr and t.t4x l0~.
=
5)< 10 ~, I .7S>( l0~
B. L. Markham, F. Rosenherger / Diffusive Pr = 0.7
0.8
‘
0.6
Sch
I
))
= 1.0
1/
-______ —
0.4
—
___________
-
/
I
—
convective vapor transport
1.0
0.8
-
/
J
‘
-
I
2~1O~
~‘)
).....s~o3
,,~.._-“
-
/ /
0.6’
5
~—
G,
247
1.78x10’
0.4
-
Gr ~1.l4x10°
0.2
-
1
1.78,10’ -~
0.2
—
I ‘
1.14,06
0
‘
1
DIMENSIONLESS VELOCITY
Fig. 6. Profiles of velocity in transport direction in plane 3, 1.78 X 1O’~and 1.14x 106. x = L/2 of cavity for Gr = 5 x tO Note that velocity is not zero at mid-height.
cantly below midheight. Also, the two lobes of the profile have different areas. This asymmetry has two reasons:
Sch
I
2
1.0
3
4
DIMENSIONLESS GROWTH FLUX
Fig. 7. Dimensionless flux into growing crystal at Gr = 2>< io~, 5 X iO~,1.78 x jØ4 and 1.14 x 106. Flux for 1,14 x 10~’has been divided by 3. Maximum in flux for Gr = 1.78 x 10~is equivalent to a local growth rate of 1 mm/h in a system in which component A is iodine. Dashed line represents Gr = 0 case. Flux made dimensionless by dividing by pAP/h, where PA is the
density of A at the interface.
(a) the net mass transfer across the cavity; more
mass is transported “down the tube” than flows back to the source;
convection on the overall transport rate is small, convection causes strong deviations from uniform
(b) the thermal expansion of the vapor; the vapor
species distribution. One must realize, however, that even mere diffusion in viscous interaction
flowing in the top part of the cavity is warmer and, hence, less dense than the vapor in the lower part. In addition to the asymmetry in the mid-length plane the whole convection role tends to tilt with respect to the side wall such that the y-location of the cross-over referred to above increases with x. These details cannot be obtained from “incompressible” descriptions (i.e. p const.) such as the
with the side walls, can result in some non-uniformity in growth flux [3]. Yet, for the underlying case this non-uniformity is only a few percent. Hence, the major deformation in the flux profile observed here results from free convection. Note that for Gr 5 x I0~the flux variation across the interface amounts to a factor of 5 and at Gr 1.78 x iO~to a factor of about 20. At Gr 1.14>< 106,
KU model,
the variation is reduced again. This, however, is a
Of particular importance to the crystal grower is the growth flux distribution across the interface, In fig. 7 we have plotted the component velocity
computational artifact due to the finite grid size used. At very high Gr’s, the major concentration
=
=
=
=
for A on the crystal interface (“growth flux”) as it results for our system for AR = 5 at four different Or values. As can be expected one finds in general a higher nutrient supply in the upper part of the cavity, where the flow connects source and crystal
change takes place over a distance comparable to our grid spacing so that the numerical results become unrealistic. This “boundary layer” development is also reflected in a shift towards the upper wall of the maximum flux region, as can be seen from fig. 7. This is associated with the tilt of
more directly. Note also that in the lowest part, the flux is reduced below its “diffusion-only” value (dashed line). The magnitude of the uniformity is surprising. Even for Or 2 x iO~,where Sh 2, the flux varies by a factor of 2.4 across the interface. Hence, even when the net effect of free
the convection roll mentioned before. These specific non-uniformities in interfacial flux were obtained for Pc 0.94 and Sch 1.0. Other combinations of these dimensionless groups will lead to different profiles, though the general tend will remain the same. For systems with lower
=
=
=
248
B. L. Markham, F. Rosenherger
Pc (lower growth rate) the magnitude of the nonuniformity will be reduced. A lower Sch (e.g. higher binary diffusivity) will have the same effect. 3.2. Effects of Prandt/ and Schmidt number Since KU indirectly assumed a Pr —s 0 throughout their analysis we have investigated the effect of Pr variations on transport. It is common practice to assign a Pr = 0.7, i.e. the value for air, to VT systems. We have estimated, though, that due to the wide variation of vapor composition used in the laboratory practice a representative range for VT is 0.4 ~ Pr ~ 1.5. To well cover this range together with KU’s idealized assumption, we have used a span of 0 ~ Pr ~ 2 in our calculations. Fig. 8 shows the results obtained for Or = 7.1 X 10~for the two aspect ratios used above. For Pr = 2 (not shown in fig. 8) and AR = 5 we found a Sh = 13.5. It is interesting to note that the maximum overall transport is not obtained at Pr = 0,
where (as a consequence of
it
—~
oc) the tempera-
ture gradient in the core remains unchanged. For AR = 5, the maximum Sh lies at Pr 0.25 and for AR = 10 at Pr 0.1. This well demonstrates that
/
Diffusive
—
convective (roper transport
the net transport between source and crystal in the convection regime results from a precarious halance between core-driven and boundary-layer (end region) driven flow. As depicted in fig. 9 for AR = 5, the longitudinal temperature profile (plotted for v = h/2) for Pr = 0.1 deviates in a large part of the core only little from the linear profile (dashed) of a mere conduction state (Pr = oc). Hence, the flow will be dominantly coredriven. At Pr = 0.3, one finds a significant temperature drop in the end regions. Hence, the core and boundary layer mechanisms will contribute with comparable strength to the overall flow. Finally, the temperature profile for Pr = 2 suggests dominant boundary layer driven flow. This is also supported by the change in the interfacial flux distributions plotted in fig. 10. In spite of the relatively low Or (compared with the highest case in fig. 7) the change from Pr = 0.1 to Pr = 2.0 causes a significant shift of the growth flux maximum towards the upper boundary. For comparison we have also plotted the KU predictions in fig. 8. Note that even at Pr = 0. that was assumed in their analysis. there is a large discrepancy to our Sh values. This illustrates again that the other simplifying assumptions made are as
important for the limitations of the KU model as 26
the constant temperature gradient.
~,,-
For a further characterization of the trends in VT flows we have also varied the Schmidt number
24
~
22
-
20
-
in our model system. The Schmidt number conAR
-
I0
~
1~: ~
12
______________
~
AR -10
0
0.2
1 I 1.1,1 0.4 0.6 0.8 1.0 PRANDTL NUMBER
1.2
: ~=:*~
0.2-
A65
I 1.4
~T’
Fig. 8. Dependence of Sherwood number on Prandtl number for AR = 5 and 10. Heavier lines: our model lighter lines: KU.
-
:
0
0
Sch~1.O
I
0.1
I 0.2
I 0.3
I
I
I
0.4
0.5 xIL
0.6
Fig. 9. Temperature profiles in plane v and 2.0.
I
=
0.7
0.8
0.9
1.0
h /2 for Pr = 0.1. 0.3
B.L. Markham, F. Rosenberger
)Sch ARI
0.8 1.0
= =
1.0I 5
I
0.6
0.4
I
I
=
‘
/
fleets the12width of the momentum concentrationand 11 boundary display layers our(at results the interface) for and the Sh [20]. depenFigs. dence on Sch for AR = 5 and 10, respectively, computed for fixed Gr = 7.1 X iO~and Pr = 0.7. The corresponding KU results are also given in
.
‘
these figures. One sees that at low Sch values, -
port. KU’s At analysis high tends Schmidt to overestimate numbers (relatively the net translow binary diffusivities) their Sh values saturate, whereas ours increase steadily with Sch.
0.2
0
,
0
1.0
I 2.0
I 3.0
4.0
5,0
6.0
3.3. Effects of geometry
DIMENSIONLESS GROWTH FLUX
Fig. 10. Dimensionless flux into growing crystal for Pr =
0.1,
3.3.1. Aspect ratio
0.3 and 2.0.
35
1
I
Gr
7,lx 113’
AR
5
I
I
I
I I
1.6
1.8
To demonstrate the effect of the aspect ratio on transport we have kept all boundary conditions and the height h constant and varied only the length of the horizontal enclosure. Fig. 13 shows the results for a system that is identical to those treated in our earlier work on vertical ampoules [10] (Pc = 0.94, Pr = 1, Sch = 0.15, h = 2 cm) for Or = 1.78 X iOn. Also plotted is KU’s prediction. The two curves a have similar shape with maxima at AR 3. There is good agreement for AR = 10. For lower AR values, where their asymptotic stream function assumption becomes unrealistic,
25
15
249
Diffusive — convective vapor transport
-
0
0.2
0.4
0.6
0.8 1.0 1.2 SCHMIDT NUMBER
1.4
2.0
Fig. 11. Dependence of Sherwood number on Schmidt number at AR = 5. Heavier line: our model, lighter line: KU.
I
1.8
I I
1.4
1.6
1.8
-
2,0
SCHMIDT NUMBER
Fig. 12. Dependence of Sherwood number on Schmidt number with AR = 10. Heavier line: our model, lighter line: KU.
wiGz 1.0 0 a 14 0 12
I 4
-
Gr=1,78x10 Pr = 1.0 Sch = 0.1 5
-
_________________________________________________________ I I I I I I I I
w Gr=7.1x10’ Pr=0.7 = ~20-30 -/111111 0 0________________________ 0.2 0.4 0.6 0.8 1.0 1.2
I
,
2I
4I 6I ASPECT RATIO
8I
10I
-
Fig. 13. Dependence of Sherwood number on AR for same conditions as in ref. [3]. Heavier line: our model, lighter line: KU.
B. L. Markham, F. Ro.cenherger
250
I
~
I
Sch
30
I
/
az
convective (roper transport
to gravity. We use the angle 4) between the main transport direction (x-axis) and the gravity vector (fig. 2) to characterize the inclination. Hence 4 =
-
90° corresponds to the horizontal systems dis-
w
0 ~ w 0 I
—
of the orientation of the VT enclosure with respect I
1.0
-
Diffusive
cussed above, and 4) =
~ 20 25 is 10-
~-~56xI~~
-
and 180°. respectively.
00
represent the Pr “top-heated” “bottom-heated” situations. ion. ForWe comparison have chosen with AR our = and earlier 3 i.e. simulations the vertical system with ampoules = 1. [101 Sch we = have 0.15. first and considered Ormaximum = 1.78for ~a
5
transport situation in fig. 13. Note that a fixed (Jr U
~
0
I
I
2
3
I
4 5 6 ASPECT RATIO
7
8
9
10
Fig. 14. variation of Sherwood number with aspect ratio at 4, 3.56x 10~and 2.85x iO~.Heavier lines: our GR =1.78x10 model, tighter lines: KU.
KU’s Sh values are consistently higher than ours.
If one changes the Schmidt number to 1.0 and the Prandtl number to 0.7 the Sh(AR) function changes significantly for the same Or (see Or = 1.78 x iO~in fig. 14). The increase in transport (Sh) scales roughly with the decrease in diffusivity, which was changed by about a factor of 6 to obtain Sch = I rather than 0.15. One must recall that the Sh is normalized with respect to the mere diffusive flux. As DAB drops so does the diffusive flux and convection becomes relatively more effective in increasing Sh. Besides this change in Sh, also the maximum of Sh(AR) shifts from 3 to
statement with variation of the orientation is physically rather meaningless and should only indicate that the vapor pressure and thermal boundary conditions remained unchanged during the changes in inclination. Fig. 15 shows the Sh dependence on 4) for the low Sch system. Results for Pr = 0.7, Sch = 1 and a more commonly used AR = 10 are given in fig. 16. It is interesting to see that maximum transport does not occur in the bottom-heated configuration, where the destabilizing temperature gradient is parallel to the gravity vector. For both systems (figs. 15 and 16) the maximum Sh lies at 120°---130°.This is in good agreement with the
~T
-~
I.s~f
-
1 --
.
Or =i.78x10’ Pr —1 ij 1.4
--
approximately 5. Fig. 14 shows also that the loca-
1
ScIs 015 AR=3
tion of the maximum sensitively depends on Or
uJ
With an increase to Or = 3.56 x IO~, the maximum shifts to AR = 7. For Or = 2.85 x iO~,i.e. at
the onset of boundary layer behavior (see, e.g., fig. 5), no maximum occurs in the range plotted. For
w
1.3
z a 0
I I
1.2
-
-~
UI
fully developed boundary layer flow one can cxpect that Sh is proportional to the ratio L over concentration boundary layer width 6~.Since ~ (Or)i/4 [9,19], no maximum should be observed. For comparison we have also plotted the corresponding KU curves (light lines) in fig. 14.
3,3.2. Inclination Now we will discuss the results from variations
1,.
1.0 -
0’
20’
40’
60’ 80’ 100’ 120’ INCLINATION ~s
140’ 160’ 180’
Fig. 15. Dependence of Sherwood number on inclination at AR 3. Note that maximum occurs between 9O0~I8(I0. i.e. not for the vertical, bottom-heated enclosure.
B. L. Markham, F. Rosenberger I
-
I
F
/
field on orientation change. In horizontal en-
~TT~T~
Pr =0.7 Sch =1.0 AR =10
more or less parallel to the side walls. On inclination, however, the convection roll becomes significantly tilted with respect to the side wall, reflecting the complex interaction of horizontal and vertical temperature (density) gradients in these systems. closures (4) = 90°) the core flow (fig. 6) proceeds we found (after considerable computational diffiFor the bottom-heated configuration (4) = 180°)
. -
6
-
8 ~ O~ 4 w I UI
culties) a complete change flow character. The asymmetric unicellular roll, in schematically depicted
3
-
0..r 0’
251
Diffusive— convective vapor transport
I 20’
I
40’
I I I J~,..I I 60’ 80’ 100’ 120’ 140’ i60’lBO’ INCLINATION
~
Fig. 16. Dependence of Sherwood number on inclination at AR = 10.
mere convective behavior in rectangular enclosures. Earlier numerical results show a maximum in the Nusselt number (thermal Sherwood number) at similar values [9,21,22]. Schinkel [9] has found for AR ~ 1, that the maximum Nusselt number is reached at 1200. For systems with adiabatic walls but AR’s and Or’s similar to ours [21,22] the Nusselt number has been found to reach its largest value at 130°.It appears that only Wirtz and Tseng [21] have considered a systern with conducting walls and an AR similar to ours. Their solutions indicate that the highest Nusselt number is reached at 1350. Suprisingly how—
—
ever, considering our work and that of Schinkel, these authors found that the Nusselt number stayed at anemometry constant studies value showed for larger also 4)’s. that Laser forVT the above Or maximum occur 900 <4) atthe < AR around transport 1800 the was maximum 4)vector = rates also 130° seen atwhich convection [14]. some in actual This angle trend velocities between experifor ments based on (ref. an [7], effective fig. 2). Or Hence, in transport simply aDoppler analysis projectiona and of gravity onto the main transport direction is used (eq. (55) in ref. [7]) cannot be expected to yield realistic results. The complexity of the Sh results are well reflected in the changes that occur in the velocity
in fig. 1, gave way to a double roll with symmetry about the middle plane y = h/2: upward flow in the middle of the enclosure and downward recirculation along the side walls. This is illustrated in fig. 17 by a U-velocity profile at x = L/2 for Or = 1.78 X 10”. Note that the area of the positive velocity exceeds those of the negative (downward)
regions. This is due to the net transport of A across (up) the cavity. For 4) = 00, i.e. for top-heating, we found also symmetry of the velocity profile with respect to the middle plane y = h/2, as depicted in fig. 17. This curve clearly shows the presence of convection even for top heating. Diffusive—advective mass transport from the hotter source tends to raise the temperature along the midplane more than closer to the (thermally stabilizing) side walls. Hence, horizontal temperature gradients and expansive convection develops, which leads to the character-
~
‘~
~
0,J _____ 0.6 0.4 0.2 0 -0.1
11” 0
I 0.1
~
180 Or Pr—0.7 =1.78x10’ Sch ~R =1.0 =10 0.2 0.3
DIMENSIONLESS VELOCITY
Fig. 17. Profile of velocity in transport direction in plane x = L/2 at 4 = 0 and 1800, showing symmetry about plane y = h/2.
252
B. L. Markham. F. Ro.senherger
/
istic dip in the velocity profile in fig. 17 (4) = 0). This we have discussed before in more detail in connection with solutal convection [10,23]. Of course, in treatments with uncoupled mass, momentum and energy transport, such insight cannot be obtained. Here we should also point out that a velocity profile similar to the 4) = 0 case of fig. 17 has been observed previously [18]. However, the authors interpreted this profile merely in terms of a diffusive deformation of a parabolic flow profile rather than the convective flow that results from the horizontal density gradients in such a configuration. The observed flow symmetry with respect to = h/2 is of some fundamental interest. For mere convection in vertical cylinders, rolls at not too high Or are asymmetric. In our earlier work on VT in vertical cylinders, however, we had assumed axisymmetric flows (for numerical simplicity’s sake) and argued that these will prevail due to the (symmetric) net transport in the system [101. The above findings give limited (i.e. 2D) support to this assumption. Finally, we have investigated the growth flux non-uniformities that arise in inclined cavities. In general, for 4) < 90° the non-uniformity is somewhat reduced as compared to fig. 7. For 4)> 90° it is first slightly increased up to the maximum transport angle, 4)max and then reduced. In fig. 18 we compare the growth flux distribution for AR = 10, Or = 1.78 X iO~and 4)fu~ = 130° (see fig. 16) with
~
1.0
\
Or =i.78xi0’ 0.8
/
AR =10 06
~/ ,,/
-
90~~>/
/~
0.4
0.2
/ 7 /
-
f 0
\ )
~ ,/
~
_•_7~
0.6
0.8
1.0
lower (upper) half of the enclosure. The remaining part changes the flow direction by dropping (rising) at the side wall over a considerably distance from the respective end. Further evidence of the fundamental shortcon3ings of 2D discussions comes from current 3D
below our model predictions.
1.21.4
DIMENSIONLESS GROWTH FLUX
Fig. 18 Dimensionless growth flux distribution for and 130°.
In order to compare our predictions with cxperimental results we have performed a detailed literature search, There are, however, only few VT experiments that have been characterized to an extent that encourages a comparison [6,7]. Even with these experiments we can only obtain semiquantitative correlation. One could list numerous parameters that are not accounted for even in our computationally complex model and that may he important in reality. Among these are the neglect of radiative heat transfer that tends to reduce the effective Or, the often more complicated interface geometries obtained in reality and the neglect of the Soret effect. Yet, we believe that the lack of quantitative agreement has more fundamental rcasons, Recent laser Doppler anemometry studies of the convective velocity field in cylinders performed in our laboratory have revealed rather unexpected 3D flow features [13]. It was found that only part of the vapor flows along the cold (hot) thermode. i.e. the “crystal” and “source” in fig. I, into the
functions of figs. 3 and 4 ought to be shifted to the right. This trend is supported by the observation that most experimental results for Sh appear to lie
,/__
0.4
3.4. Comparison with experimenl.v
of the Fluid Dynamics Oroup at the University of Marseille and our group [15,16]. These efforts show that, in general, 2D systems enter the various fluid dynamics regimes. (such as boundary layer flow) at lower Or values than the corresponding 3D systems. This implies, for instance, that the Sh(Gr) .
/
02
the distribution obtained for 4) = 90°. Besides the increase in transport one notes also a shift in the location of the flux maximum as a result of the before mentioned tilt of the convection roll.
simulation that is being performed in a joint effort
\ \
-
Diffusive — conI’ectil’e vapor tra/isporl
~
=
900
Why then strive for improved 2D models, as we have attempted in this work? Full-fledged 3D treatments of the three coupled transport equa-
B.L. Markham. F. Rosenberger
/
tions for mass, momentum and energy (including
concentration-induced density gradients) for realistic boundary conditions will remain beyond reach for most workers for some time to come. Hence, 2D models are currently the only theoretical means of obtaining insight on practically important trends that can lead to improvements of the transport conditions.
4. Summary Our numerical solutions for horizontal VT agree well with predictions of the KU model at low Grashof numbers and large aspect ratios. At medium and high tend to level off. Or, Our KU’s modelSherwood yields a numbers realistic increase which asymptotically approaches a (Gr)°25 dependence, in agreement with boundary layer theory. At low aspect ratios and moderate products of Sch and Or, KU’s analysis tends to overestimate the Sh. In both approaches a maxi-
Diffusive— convective vapor transport
NSO-1534 and by the National Science Foundation under grant DMR-79-13183 are gratefully
acknowledged. Stimulating discussions with P. Bontoux are also very much appreciated.
Nomenclature A AR
Gaseous growth component Aspect ratio (L/h)
B C 1,
Gaseous inert component
DAB g Or h k L p 1 P Pc
mum in Sh is reached at smaller aspect ratios, the exact value of AR slowly growing with Or as long as full boundary layer behavior is not achieved. The Schmidt and Prandtl numbers of the vapor were shown to have significant influence on transport, effects that cannot realistically be deduced from the simplifying KU model. We have also investigated the effect of orientation of the enclosure on Sh. Transport becomes maximum for an angle of about 300_400 between the main transport direction and the horizontal with the hotter end below. The exact value of this angle, however, must be expected to depend on Or and Ar. Beyond the overall transport predictions, our analysis yields also solutions for temperature and growth flux distributions. These show practically important non-uniformities, and shed light on the physical origin of shortcomings of earlier transport models.
Specific heat at constant pressure per unit mass Binary diffusion coefficient Acceleration due to gravity 3~T/~2) Orashof Height ofnumber cavity ($gh Thermal conductivity Length of cavity Partial pressure of ith (A or B) component Total pressure Péclet number (ln[ ~vB( L )/p ~~(0)], see ref.
[3]) Pr Sch Sh
T L~T T
Prandtl number (v/K) Schmidt number (v/DAB) Sherwood number (ratio of total interfacial transport fluxes in 2D and corresponding 1D model) Temperature T~— T~
Normalized temperature [(T T )/( T~ 1’)] —
7~
T~ u
U v V
W~ x X y
Acknowledgments Y Support of this work by the National Aeronautics and Space Administration under grant
253
—
Temperature at crystal Temperature at source Component of velocity in x-direction Dimensionless x-component of velocity (uh/v) Component of velocity in y-direction Dimensionless y-component of velocity (vh/v) Weight fraction of I th (A or B) component Coordinate in main transport direction Dimensionless x-coordinate (x/L) Coordinate perpendicular to main transport direction Dimensionless y-coordinate (y/h)
254
1k!.. Markham. F Rosenherger
/
DiJfuso’e-— ionIeelIvr’ vapor trailsport
Greek letters
Enclosures (Dutch
Efficiency Bureau-- Pijnacker.
t)elft.
1980).
,8
Thermal expansion coefficient (ideal gas:$ = 1 /T / . . . Thermal diffusivity Viscosity Kinematic viscosity Mass density
t~
p. v
p
12]
Crystal Growth 51(1981) 426. 111 R. Jackson. Transport iii Porous Amsterdam, 1977) eli. 4.
Catalysts
Elsesier,
112) B.S. Jhaveri. B.L. Markham and F. Rosenherger. (heni. Eng. Comntun. 13(1981) 65.
[131(LII.
Schirok~and F. Rosenherger. Intern. .1. Heat Mass Fransfer 27 (1984) 587, [14] F. Rosenberger and AC. 1 Iurford. in preparation. 115] C. Smutek. P. Bontoux, B. Rou\. (iii. Schiroky, AC. Hurford. F. Rosenherger and U. de VahI Davis, Nunieri-
References )l] H. Wiedemeier. F.C. Klaessig. E.A, Irene and Si. Wev. Crystal Growth 31 (1975) 36.
101 BE. Markham, [).W. Greenwell and F. Rosenhcrgcr, J.
j.
H. Wiedemeier. J. Sadeek, F.C. Klaessig, M. Norek and R. Santandrea, J. Electrochem. Soc. 124 (1977) 1095.
cal Heat Transfer, submitted.
116] P. Bontoux. B. Roux, G.H. Schiroky, B.L. Markham and F. Rosenberger, Intern. J. 1-leat Mass Transfer. submitted. )17[
A.P. Gosman and W.M. Pun. Report No. UTS/74/2. Impenal College of Science and Technology, London,
[3] D.W. Greenwell, B.L. Markham and F. Rosenherger. J. Crystal Growth 51(1981) 413. [4] A. Solan and S. Ostraeh, Convection Effects in Crystal Growth by Closed-Tube Chemical Vapor Transport, in: Preparation and Properties of Solid State Materials. Vol.
1973. [181 G. Wyatt. W.E. Ihele and ER. Eckert. Intern. .1. Heat Mass Transfer 14(1973) 513. [19) A.E. Gill, J. Fluid Mech. 26 (1966) 515. [20] F. Rosenberger, Fundamentals of Crystal Growth I
4, Ed. W.R. Wilcox (Dekker, New York. 1979) pp. 63--I 10. [5] K. Klosse and P. Vllersnia, J. Crystal Growth 18 (1973) 167. [6) C’. Paorici and C. Pelosi. J. (‘rystal Growth 35(1976) 65.
(Springer, New York. 1979) p. 300. 1211 R.A. Wirtz and W.F. Tseng, Natural Convection Across Tilted Rectangular Enclosures of Small Aspect Ratio, in: Natural Convection in Fnclosures—EITD. Vol. 8. Eds. J.
[71 D.
Chandra and H. Wiedemeier, J. Crystal Growth 57 (1982) 159. [8] B.S. Jhaveri arid F. Rosenberger. J. Crystal Growth 57
Catton and K.E. Torrance (ASME, New York. 198W pp. 47—54. 22] iC. Grondin, Phi) Thesis. Université de Prosence.
(1982) 57. [9] W. Schinkel, Natural Convection in Inclined Air-Filled
Marseille (1978). [23) B.L. Markham and F. Rosenberger, Chem. Eng. (‘ornniun. 5(1980) 287.
241
DIFFUSIVE-CONVECTIVE VAPOR TRANSPORT ACROSS HORIZONTAL AND INCLINED RECTANGULAR ENCLOSURES B.L. MARKHAM and F. ROSENBERGER Department of Physics, University of Utah, Salt Lake City, Utah 84112, USA
Received 23 August 1983; manuscript received in final form 13 February 1984
The enhancement of vapor transport across horizontal and inclined rectangular (2D) enclosures by expansive convection is studied numerically for a range in Grashof number that extends up to boundary layer flow. Comparison is made with predictions of the simplifying model and resulting analytical treatment of Klosse and Ullersma for various aspect ratios of the enclosure. The influence of varying Schmidt and Prandtl numbers is investigated. Solutions for interfacial (growth) flux and temperature distributions are given, showing practically important non-uniformities.
I. Introduction There has been considerable interest in accurate predictions of transport rates for crystal growth from vapors in closed ampoules. In particular, this interest has grown since chemical vapor transport experiments recently conducted under low-gravity conditions aboard space-crafts gave results [1,2] that still defy quantitative interpretation. The earlier, mostly 1D or idealized 2D treatments of vapor transport (VT) in closed ampoules have been reviewed elsewhere [3,4]. Among these, the most realistic and, hence, widely used model for diffusive—convective horizontal VT is due to Kiosse and Ullersma [5] (hereafter referred to as KU). Various experimental studies, however, failed to quantitatively confirm KU’s predictions for a significant range of characteristic parameters; see, e.g., refs. [6,7]. Most of the simplifying assumptions made by KU to obtain a closed form analyti cal solution were relaxed in later numerical work [8]. Due to numerical difficulties, however, the Boussinesq approximation was used and only flows with rather low Grashof number were considered. For this limited parameter range good agreement with the simpler KU model was obtained [8]. In the present paper, based on an improved numerical approach, we cover the whole convection
(Grashof number) range that is of interest in the laboratory. In addition to the horizontal transport problem we consider the inclination of the rectangular enclosure with respect to gravity over the whole range between “top-heating” and “bottomheating”. Stimulated by some of the simplifications in KU’s work, we have also explored the influence on transport rates of variations in vapor properties (Schmidt and Prandtl numbers) over a range that is representative of vapor mixtures in the VT practice. Since KU’s work forms the starting point for our efforts we shall give a brief summary of their model and approach for solving it. A rectangular horizontal cavity of height h and horizontal extent L is assumed, see fig. 1. The third dimension is taken as very large. The vertical end walls are ___________
L
-___________________________
h T5
T C
Fig. 1. Expansive convection due to horizontal temperature gradient in rectangular enclosure. Overall flow, and velocity profile (S-profile) in vertical Center plane; schematic.
0022-0248/84/$03.OO © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
242
B, 1.. Markham. f Rosenherger / l)s,/,fuses ‘e --~‘ons’ee1,s’elapor transport
differentially heated to 1~and J~.and drive expansive convection of the vapors between them. The end walls form also the interfaces between solid source and vapor (at L,), and vapor and crystal (at 7~.),respectively. In solving for the mass transport rate, KU disregard the heat transfer across the enclosure. A linear temperature gradient aT/ax = (7’~— 7’~)/Lis assumed in the vapor, irrespective of the actual convection conditions. Note that with this assumption it is implied that the vapor has infinite thermal diffusivity and, hence, vanishing Prandtl number, Pr (see nomenclature at end of paper). Solutal density gradients are ignored. Hence, the fixed temperature gradient results also in a fixed density gradient ap/ax. Furthermore, velocity contributions arising from the net mass transfer (vapor in-flux and out-flux at source and crystal, respectively) are also ignored. The convective flow is then obtained by KU by assuming a stream function for an inifitely long, horizontal channel with anti-parallel flow (see the S-profile in fig. 1). The boundary conditions for the vertical end walls were grafted onto the stream function and the resulting unknown parameters were solved for with a variational technique. This newly found velocity distribution was then used to solve for the concentration profile. The integral solution allowed KU to make predictions about the enhancement of mass transport due to convection in the form of a rather simple mathematical function of aspect ratio AR = L/h. Grashof number Gr and Schmidt number Sch. No solutions for local concentration distribution and growth flux were given, Intuitively one can expect the following limitations from these model assumptions and the method of solution. Since a specific symmetric stream function was chosen, the model cannot account for deformations of the velocity field (development of boundary layer behavior) that arise in reality at higher flow velocities i e higher Or Synonymous with this restriction is the assumption of a constant density gradient (equivalent to Pr = 0). This restricts the drive for convection to the core mechanism. In reality, T(x) deviates increasingly from the assumed linear profile with increasing convection strength, and eventually boundary layer entrainment provides the dominant drive at
high Or (9J. For vapors, where Pr is of order one. the restriction to the core mechanism can he cxpected to he a serious limitation even at not too strong convection conditions. Also, since the stream function assumed for the core has the character of the flow in an infinitely long tube, one can expect poor results even at modest Or for systems of small aspect ratio. In the following we will quantitatively explore these questions. together with an expansion of the model to arbitrary orientations with respect to gravity, and present solutions for concentration, temperature and growth flux distributions.
2. Model and theoretical analysis We shall consider a rectangular cavity of height h and length L, the end walls of which are inclined with angle ~ + s~/2 with respect to the gravity vector, see fig. 2. The cavity is assumed sufficiently deep so that the third dimension may he neglected. As in the KU model (fig. 1) transport takes place between two flat, isothermal interfaces that are kept at 7’~,and 7., respectively, with 1~>i~.(‘omponent A sublirnes from the source, is transported through and with a non-reactive component B, and condenses at the crystal interface. We assume complete rejection of B and equilibrium of A at both interfaces. The bounding walls are assumed to be impermeable to both components and thermally ideally conducting. Radiative heat transfer is neglected. The molecular weights of the two vapor components are assumed equal, MA = M 15, which restricts our model to expansive (“thermal”) convection. Also, frictional contributions to energy
~ Source~
1
I
j
——_____
Vapors
_______
~~taI
~
~B [x
/ Fig. 2 Model for diffusive—convective vapor transport across rectangular enclosure of arbitrary orientation ~ with respect to gravity.
B. L. Markham, F. Rosenberger / Diffusive — convective vapor transport
transport are ignored. The viscosity, binary diffusivity and thermal diffusivity are assumed concentration and temperature independent, The transport in this cavity is governed by the system of nonlinear, coupled conservation equations for momentum (Navier—Stokes), mass (continuity), species (diffusion) and energy. These equations are, respectively, in rectangular coordinates and steady state
a
a
~(puu)
=
+~(pvu)
a~
F
a
au
2 /
au
av il
a2v + ~ a2u ayax —~+pg~,
ay
ay
—
rigidity v(x,0)
=
v(x,h)
=
0,
(6)
no-slip u(x,0)
=
u(x,h)
=
0,
(7)
conducting
(Ia)
a~
av
2 /
au
av
a
(lb)
0,
(2)
—
T,,—T~ L
(8)
a
au’ç
1
—
—
DAB awA(o,Y) W,~(0,y) ax
‘
(9)
v(0,y)=0,
(10)
uniform B-concentration (11)
uniform temperature
a
awç \
/
~
ay a(pvT)
=
pB(0,y)=const.,
-~—(puW,)+-~—(pvWç)
~(puT)+ a
—
—
no-slip
a2v ax
=
u(0,y)
~
+j.t—+j.t——~+pg5’,
a(pu) a(pv) ax + ay
aT(x,h) ax
—
out-diffusion of A and impermeability to B
a F
a2u axay
aT(x,o) ax Source:
a(puv) a(puv) ax + ay =
the transport region. The boundary conditions are: Walls: impermeability awA(X,o) — 8J’VA(x,h) = 0, (5)
ax +1L~
243
2T
~
k (a
a2T
T(0,y)
=
T 5.
(12)
(3)
Crystal:
(4)
in-diffusion of A and to B DABimpermeability aw~(L,y) =
-
u(L,y) As in our work on diffusive—convective transport in vertical, cylindrical ampoules [3,10], we have
kept the “compressible terms” in the momentum transport equations. Also, neither creeping flow nor boundary layer assumptions were made, in order to retain applicability to the wide range in flow parameters (Gr) to be investigated. In addition to (1)—(4), the ideal gas law and Dalton’s law of partial pressures were used to link density to temperature and concentration. To solve these equations we must state boundary conditions on both walls and interfaces enclosing
no-slip v(L,y)
1— WA(L,y)
=
0,
ax
(13)
(14)
uniform B-concentration pB(L,y) = const.,
(is)
uniform temperature T( L y) T~.
(16)
,
=
Clearly, some of these boundary conditions merit more discussion. Perfectly conducting walls
244
B.L. Markham, F Rasenherger
(or, for that matter, perfectly insulating walls) cannot be obtained in a crystal growth process. Hence, (8) must be understood as one approximation which is probably reasonable, comparing, for instance, the thermal conductivities of quartz and typical gases. The no-slip condition (7), has been
/
Diffusive
-
conteclIse taper transport
molecular weight of the vapor components any diffusion slip can safely be excluded [111. With MA = MB we can also neglect solutal convection effects. As we have pointed out before [3,12] there is a singularity at each corner of the cavity. The singularity arises from the simultaneous statement of the no-slip condition and of in- and out-diffusion, respectively, at these points. For a thorough discussion of this problem the reader is directed to ref. [12]. Suffice it to say that these singularities
out to he highly three-dimensional, even in horizontal cylinders with ideal linear temperature profile. ‘This forms probably the strongest limitation of the following results and discussion. However, the tools (codes and computer power) for a full-fledged 3D treatment of this problem including momentum, energy and species transport. are Just becoming available. Hence, the semi-quantitative guidance that can he obtained from our resuits, can probably he useful for the practitioner until more rigorous treatments can be carried out for the high Or conditions often found in the laboratory. The system of non-linear, coupled partial differential equations (1)--(4) was solved with a finite difference method. We developed an algorithm which is a variation of the code TEACH [171. The grid system used had 16 and 32 grid lines in the r-
appear to have little effect on the solutions given here. In reality the interface shape is a result of the A-flux distribution, i.e. for non-uniform interfacial
and x-direction, respectively. For some of the soiutions, in particular at higher Or, we used non-uniform grid spacings in the x-direction. with tighter grid packing close to the interfaces, where higher
flux, the shape of the interface will change and not remain flat as we have assumed. The evolution of the interface shape as a function of the prevailing
gradients in velocity, concentration and temperature are encountered. As starting solutions for the iteration process we used typically solutions to the corresponding I D problem with linear temperature gradient. However, for high Or conditions, this approach offers no advantage over an “initial fill” with zeroes and a “boot-strapping” procedure from lower to higher Gr becomes very important for efficient computation. During a computational run the total pressure was held constant only at one point, namely at the intersection of the third x- and y-grid lines (safely spaced from the corner singularity). The pressure at all other grid points was up-dated with each iteration to satisfy continuity. The criteria for convergence were an unchanging velocity field and balance of the mass flux for all cross-sections x = const. of the cavity. The calculations were carried out on a VAX 750 computer. Double precision was used to overcome certain truncation instabilities. The results presented here, including extensive program development runs, required approximately 300 CPU hours. The establishment of the physical integrity of any numerical scheme is a non-trivial task. This is particularly true for the systems considered here due to a lack of non-intrusive measurements of
the subject of some discussion, as summarized in ref. [3]. However, with the assumption of equal
heat and mass transfer conditions (which, in turn, change with the interface shape) forms a formidable mathematical task (Stefan problem). However, the rate at which the interface moves is typically more than three orders of magnitude slower than representative vapor flow velocities. Hence, the assumption of an unchanging shape of source and crystal is reasonable. To assume planar interfaces at the respective ends of the ampoule, however, is not necessarily representative of the crystal growth practice. Most workers obtain many little crystals in these ampoule regions, unless specific provisions (seeding, compacted source, etc.) for planar interfaces are made. Yet, since the purpose of this paper consists in pointing out useful trends, we feel comfortable with this mathematically greatly simplifying assumption. During the course of this work, due to input from our parallel experimental efforts [13,14] and from 3D calculations [15,16] it has become obvious that, in general, 2D treatments of VT processes in cylindrical ampoules must be rather limited. In particular at higher Gr, the convective flow turns
B. L. Markham, F. Rosenberger
/
temperature and concentration distributions in closed vapor transport ampoules. Velocity distributions have become available only recently
40
of the experimentally obtained convection behavior. These efforts, together with the results of other numerical schemes, will be reported elsewhere [16]. For the concentration distribution obtained with our scheme we gained confidence from other runs that we tailored to simulate earlier diffusive—convective experiments with gases in semi-open systems [18]. This shall be discussed in connection with fig. 17. In the following we will present the results grouped according to effects of (a) changing temperature differences and pressure, i.e. changing Grashof number; (b) changing materials properties of the vapor, i.e. changing Prandtl and Schmidt numbers; (c) changing geometry, i.e. changing ampoule length, and orientation with respect to gravity.
‘‘‘i
‘‘‘‘‘~‘‘“
Pr
[13,14]. Hence, ir order to test our scheme, we have conducted several numerical runs with single component fills, i.e. with component B only. Good
agreement was obtained with certain 2D features
245
Diffusive — convective vapor transport
30
8
= 5
Pe
~
i.56
20
0 7
1
Sch AR
-
~
—
-
•
-
-
~ ~ -
-
KU 1
-
‘‘‘‘‘13
‘‘‘‘~~~
10
I ~‘~‘‘‘‘i’~~
I
-~
106
GRASHOF NUMBER
Fig. 3. Variation of Sherwood number with Grashof number for several Pe’s at AR = 5. Heavy line from our model calculation, lighter line prediction from KU model.
Fig. 3 illustrates, for AR = 5, the importance of convection for mass transport in VT. For Or ~ 2 x i0~, diffusion dominates, as reflected by
3. Results and discussion
60
The following dimensionless groups are stated with each result to define the system: Grashof number Gr, Prandtl number Pr, Schmidt number
50
Sch, Péclet number Pc, aspect ratio AR. In addition the angle 4) (see fig. 2) is given. The density, that enters Or, Pr and Sch, was taken at the pressure monitoring point for the average temperature(1,+7~)/2.
40
-
-
30
-
-
3.1. Effects of Grashof number
20
-
‘I
-
‘‘‘‘I
‘‘‘‘‘I’’’!’
~
-
AR
=
KU
These solutions are all for “horizontal ampoules”, i.e. for 4) 90°.Aspect ratios of 5 and =
10 are considered. Constant values for Pr
=
and Sch 1.0 were used throughout this section. Most of the Gr range was scanned by changing the (monitor point) pressure. A few low Gr results were obtained by increasing the viscosity. For comparison, a representative experimental range in Gr is from below i04—2 X 10~’[5,6]. =
10
-
-
0.7 0
.1
io~
,,,,
I
106
GRASHOF NUMBER Fig 4. Variation of Sherwood number with Grashof number for Pe 0.95 at AR 10. Heavy line: our model, lighter line: KU.
246
B. L. Markham, F Rosenherger / Diffusive
Sherwood numbers, Sh close to unity. For Or> 2 >< iO~, however, the convective contribution be-
comes dominant. For more slender “ampoules”, e.g. for AR = 10, as shown in fig. 4, the cross-over from the diffusion- to the convection-dominated regime occurs at a somewhat higher Gr. This is intuititively obvious, since with increasing AR there is viscous interaction of the flow with a longer sidewall, The corresponding KU solutions are also plotted in figs. 3 and 4. One sees that at Gr’s that are in the middle of the experimental range, KU’s predictions level out, whereas our less restricted model yields a continuous increase of Sh with Or, in good agreement with experimental observations [7]. The steep section of our curves (Or ~ 5 >< 10~) represent functions Sh ~ Gr°, where the power n is 0.28 and 0.25 for figs. 3 and 4, respectively. This compares very well with heat transfer results obtamed from boundary layer behavior [9,19]. We have also investigated the influence of the Péclet number on the Sh, i.e. the influence of the in-flow, and out-flow velocities from net mass transfer. A derivation of the relation for the Pc used here (see Nomenclature) is given in the appendix of ref. [3]. Somewhat surprising to us, a variation of 0.95 ~ Pc ~ 0.56 (which well covers the experimental range) has practically no result on Sh (fig. 3). This reflects the negligible effect of the diffusive through-flow on total velocity even in the presence of only weak convection, The corresponding solutions for the velocity field show that even at Gr = I0~,where convection aids net transport very little (i.e. Sh 1), the convective velocities are two orders of magnitude larger than the diffusive velocities. Because of this insensitivity of the results to reasonable Pc variations we have used a constant Pc = 0.95 for all following caiculaAt this point one may be tempted to neglect the species diffusion equation altogether during solving for velocity and temperature distributions and, then to determine the concentration distribution based on these uncoupled solutions. This would reduce the computation time significantly. One could even get fairly reliable information on the local fluxes at the interfaces in this way. However,
this procedure neglects whatever effect these in-
convective vapor transport
—
and out-fluxes might have on the overall flow. As we shall see, the interfacial fluxes can change the entire form of the flow. Hence, the species diffusion equation must remain coupled during the solution process. Let us now concentrate on the physical reasons for the unrealistic results of the KU approach. Fig. 5 shows computed temperature profiles taken at half height (v = h/2) along the length of the cayity, for AR = 5 and three different Or’s. At Or = 5 x 1O~,where we find excellent agreement with the Sh predicted by KU, the temperature profile is practically linear. At Or = 1.78 x iO~,where the transport results of the two models begin to significantly diverge, the temperature profile begins to he deformed by convection. At Or = 1,14 X l0~, representing the highest valtie in fig. 3, and the largest disparity with KU’s prediction. the temperature profile has assumed boundary-layer type shape. Such strong deviations from the linear distribution must, of course, lead to failure of the KU model. Whereas in reality (and in our model) a transition from core-driven to boundary-layerdriven flow occurs, the KU model is by assumption limited to core-driven flow at all Or’s. Additional insight can be obtained from computed velocity profiles. Fig. 6 represents core flow velocities u(t’) at x = L/2 for the Or values of fig. 5. We see that only for the lowest Or one obtains a symmetric “S-profile”. At higher Or the cross-over from positive to negative velocity occurs signifi-
T
F
‘~0.8
~ ~
o.~
Or
~
.14~ -~
-
~ Pr ~
0.2
0
0
0 7
SeS
=-
1.0
AR
=
5 0.2
I
I 0.4
0.6
0.8
1.0
x/L
Fig. 5. Temperature profiles at h /2 for Gr and t.t4x l0~.
=
5)< 10 ~, I .7S>( l0~
B. L. Markham, F. Rosenherger / Diffusive Pr = 0.7
0.8
‘
0.6
Sch
I
))
= 1.0
1/
-______ —
0.4
—
___________
-
/
I
—
convective vapor transport
1.0
0.8
-
/
J
‘
-
I
2~1O~
~‘)
).....s~o3
,,~.._-“
-
/ /
0.6’
5
~—
G,
247
1.78x10’
0.4
-
Gr ~1.l4x10°
0.2
-
1
1.78,10’ -~
0.2
—
I ‘
1.14,06
0
‘
1
DIMENSIONLESS VELOCITY
Fig. 6. Profiles of velocity in transport direction in plane 3, 1.78 X 1O’~and 1.14x 106. x = L/2 of cavity for Gr = 5 x tO Note that velocity is not zero at mid-height.
cantly below midheight. Also, the two lobes of the profile have different areas. This asymmetry has two reasons:
Sch
I
2
1.0
3
4
DIMENSIONLESS GROWTH FLUX
Fig. 7. Dimensionless flux into growing crystal at Gr = 2>< io~, 5 X iO~,1.78 x jØ4 and 1.14 x 106. Flux for 1,14 x 10~’has been divided by 3. Maximum in flux for Gr = 1.78 x 10~is equivalent to a local growth rate of 1 mm/h in a system in which component A is iodine. Dashed line represents Gr = 0 case. Flux made dimensionless by dividing by pAP/h, where PA is the
density of A at the interface.
(a) the net mass transfer across the cavity; more
mass is transported “down the tube” than flows back to the source;
convection on the overall transport rate is small, convection causes strong deviations from uniform
(b) the thermal expansion of the vapor; the vapor
species distribution. One must realize, however, that even mere diffusion in viscous interaction
flowing in the top part of the cavity is warmer and, hence, less dense than the vapor in the lower part. In addition to the asymmetry in the mid-length plane the whole convection role tends to tilt with respect to the side wall such that the y-location of the cross-over referred to above increases with x. These details cannot be obtained from “incompressible” descriptions (i.e. p const.) such as the
with the side walls, can result in some non-uniformity in growth flux [3]. Yet, for the underlying case this non-uniformity is only a few percent. Hence, the major deformation in the flux profile observed here results from free convection. Note that for Gr 5 x I0~the flux variation across the interface amounts to a factor of 5 and at Gr 1.78 x iO~to a factor of about 20. At Gr 1.14>< 106,
KU model,
the variation is reduced again. This, however, is a
Of particular importance to the crystal grower is the growth flux distribution across the interface, In fig. 7 we have plotted the component velocity
computational artifact due to the finite grid size used. At very high Gr’s, the major concentration
=
=
=
=
for A on the crystal interface (“growth flux”) as it results for our system for AR = 5 at four different Or values. As can be expected one finds in general a higher nutrient supply in the upper part of the cavity, where the flow connects source and crystal
change takes place over a distance comparable to our grid spacing so that the numerical results become unrealistic. This “boundary layer” development is also reflected in a shift towards the upper wall of the maximum flux region, as can be seen from fig. 7. This is associated with the tilt of
more directly. Note also that in the lowest part, the flux is reduced below its “diffusion-only” value (dashed line). The magnitude of the uniformity is surprising. Even for Or 2 x iO~,where Sh 2, the flux varies by a factor of 2.4 across the interface. Hence, even when the net effect of free
the convection roll mentioned before. These specific non-uniformities in interfacial flux were obtained for Pc 0.94 and Sch 1.0. Other combinations of these dimensionless groups will lead to different profiles, though the general tend will remain the same. For systems with lower
=
=
=
248
B. L. Markham, F. Rosenherger
Pc (lower growth rate) the magnitude of the nonuniformity will be reduced. A lower Sch (e.g. higher binary diffusivity) will have the same effect. 3.2. Effects of Prandt/ and Schmidt number Since KU indirectly assumed a Pr —s 0 throughout their analysis we have investigated the effect of Pr variations on transport. It is common practice to assign a Pr = 0.7, i.e. the value for air, to VT systems. We have estimated, though, that due to the wide variation of vapor composition used in the laboratory practice a representative range for VT is 0.4 ~ Pr ~ 1.5. To well cover this range together with KU’s idealized assumption, we have used a span of 0 ~ Pr ~ 2 in our calculations. Fig. 8 shows the results obtained for Or = 7.1 X 10~for the two aspect ratios used above. For Pr = 2 (not shown in fig. 8) and AR = 5 we found a Sh = 13.5. It is interesting to note that the maximum overall transport is not obtained at Pr = 0,
where (as a consequence of
it
—~
oc) the tempera-
ture gradient in the core remains unchanged. For AR = 5, the maximum Sh lies at Pr 0.25 and for AR = 10 at Pr 0.1. This well demonstrates that
/
Diffusive
—
convective (roper transport
the net transport between source and crystal in the convection regime results from a precarious halance between core-driven and boundary-layer (end region) driven flow. As depicted in fig. 9 for AR = 5, the longitudinal temperature profile (plotted for v = h/2) for Pr = 0.1 deviates in a large part of the core only little from the linear profile (dashed) of a mere conduction state (Pr = oc). Hence, the flow will be dominantly coredriven. At Pr = 0.3, one finds a significant temperature drop in the end regions. Hence, the core and boundary layer mechanisms will contribute with comparable strength to the overall flow. Finally, the temperature profile for Pr = 2 suggests dominant boundary layer driven flow. This is also supported by the change in the interfacial flux distributions plotted in fig. 10. In spite of the relatively low Or (compared with the highest case in fig. 7) the change from Pr = 0.1 to Pr = 2.0 causes a significant shift of the growth flux maximum towards the upper boundary. For comparison we have also plotted the KU predictions in fig. 8. Note that even at Pr = 0. that was assumed in their analysis. there is a large discrepancy to our Sh values. This illustrates again that the other simplifying assumptions made are as
important for the limitations of the KU model as 26
the constant temperature gradient.
~,,-
For a further characterization of the trends in VT flows we have also varied the Schmidt number
24
~
22
-
20
-
in our model system. The Schmidt number conAR
-
I0
~
1~: ~
12
______________
~
AR -10
0
0.2
1 I 1.1,1 0.4 0.6 0.8 1.0 PRANDTL NUMBER
1.2
: ~=:*~
0.2-
A65
I 1.4
~T’
Fig. 8. Dependence of Sherwood number on Prandtl number for AR = 5 and 10. Heavier lines: our model lighter lines: KU.
-
:
0
0
Sch~1.O
I
0.1
I 0.2
I 0.3
I
I
I
0.4
0.5 xIL
0.6
Fig. 9. Temperature profiles in plane v and 2.0.
I
=
0.7
0.8
0.9
1.0
h /2 for Pr = 0.1. 0.3
B.L. Markham, F. Rosenberger
)Sch ARI
0.8 1.0
= =
1.0I 5
I
0.6
0.4
I
I
=
‘
/
fleets the12width of the momentum concentrationand 11 boundary display layers our(at results the interface) for and the Sh [20]. depenFigs. dence on Sch for AR = 5 and 10, respectively, computed for fixed Gr = 7.1 X iO~and Pr = 0.7. The corresponding KU results are also given in
.
‘
these figures. One sees that at low Sch values, -
port. KU’s At analysis high tends Schmidt to overestimate numbers (relatively the net translow binary diffusivities) their Sh values saturate, whereas ours increase steadily with Sch.
0.2
0
,
0
1.0
I 2.0
I 3.0
4.0
5,0
6.0
3.3. Effects of geometry
DIMENSIONLESS GROWTH FLUX
Fig. 10. Dimensionless flux into growing crystal for Pr =
0.1,
3.3.1. Aspect ratio
0.3 and 2.0.
35
1
I
Gr
7,lx 113’
AR
5
I
I
I
I I
1.6
1.8
To demonstrate the effect of the aspect ratio on transport we have kept all boundary conditions and the height h constant and varied only the length of the horizontal enclosure. Fig. 13 shows the results for a system that is identical to those treated in our earlier work on vertical ampoules [10] (Pc = 0.94, Pr = 1, Sch = 0.15, h = 2 cm) for Or = 1.78 X iOn. Also plotted is KU’s prediction. The two curves a have similar shape with maxima at AR 3. There is good agreement for AR = 10. For lower AR values, where their asymptotic stream function assumption becomes unrealistic,
25
15
249
Diffusive — convective vapor transport
-
0
0.2
0.4
0.6
0.8 1.0 1.2 SCHMIDT NUMBER
1.4
2.0
Fig. 11. Dependence of Sherwood number on Schmidt number at AR = 5. Heavier line: our model, lighter line: KU.
I
1.8
I I
1.4
1.6
1.8
-
2,0
SCHMIDT NUMBER
Fig. 12. Dependence of Sherwood number on Schmidt number with AR = 10. Heavier line: our model, lighter line: KU.
wiGz 1.0 0 a 14 0 12
I 4
-
Gr=1,78x10 Pr = 1.0 Sch = 0.1 5
-
_________________________________________________________ I I I I I I I I
w Gr=7.1x10’ Pr=0.7 = ~20-30 -/111111 0 0________________________ 0.2 0.4 0.6 0.8 1.0 1.2
I
,
2I
4I 6I ASPECT RATIO
8I
10I
-
Fig. 13. Dependence of Sherwood number on AR for same conditions as in ref. [3]. Heavier line: our model, lighter line: KU.
B. L. Markham, F. Ro.cenherger
250
I
~
I
Sch
30
I
/
az
convective (roper transport
to gravity. We use the angle 4) between the main transport direction (x-axis) and the gravity vector (fig. 2) to characterize the inclination. Hence 4 =
-
90° corresponds to the horizontal systems dis-
w
0 ~ w 0 I
—
of the orientation of the VT enclosure with respect I
1.0
-
Diffusive
cussed above, and 4) =
~ 20 25 is 10-
~-~56xI~~
-
and 180°. respectively.
00
represent the Pr “top-heated” “bottom-heated” situations. ion. ForWe comparison have chosen with AR our = and earlier 3 i.e. simulations the vertical system with ampoules = 1. [101 Sch we = have 0.15. first and considered Ormaximum = 1.78for ~a
5
transport situation in fig. 13. Note that a fixed (Jr U
~
0
I
I
2
3
I
4 5 6 ASPECT RATIO
7
8
9
10
Fig. 14. variation of Sherwood number with aspect ratio at 4, 3.56x 10~and 2.85x iO~.Heavier lines: our GR =1.78x10 model, tighter lines: KU.
KU’s Sh values are consistently higher than ours.
If one changes the Schmidt number to 1.0 and the Prandtl number to 0.7 the Sh(AR) function changes significantly for the same Or (see Or = 1.78 x iO~in fig. 14). The increase in transport (Sh) scales roughly with the decrease in diffusivity, which was changed by about a factor of 6 to obtain Sch = I rather than 0.15. One must recall that the Sh is normalized with respect to the mere diffusive flux. As DAB drops so does the diffusive flux and convection becomes relatively more effective in increasing Sh. Besides this change in Sh, also the maximum of Sh(AR) shifts from 3 to
statement with variation of the orientation is physically rather meaningless and should only indicate that the vapor pressure and thermal boundary conditions remained unchanged during the changes in inclination. Fig. 15 shows the Sh dependence on 4) for the low Sch system. Results for Pr = 0.7, Sch = 1 and a more commonly used AR = 10 are given in fig. 16. It is interesting to see that maximum transport does not occur in the bottom-heated configuration, where the destabilizing temperature gradient is parallel to the gravity vector. For both systems (figs. 15 and 16) the maximum Sh lies at 120°---130°.This is in good agreement with the
~T
-~
I.s~f
-
1 --
.
Or =i.78x10’ Pr —1 ij 1.4
--
approximately 5. Fig. 14 shows also that the loca-
1
ScIs 015 AR=3
tion of the maximum sensitively depends on Or
uJ
With an increase to Or = 3.56 x IO~, the maximum shifts to AR = 7. For Or = 2.85 x iO~,i.e. at
the onset of boundary layer behavior (see, e.g., fig. 5), no maximum occurs in the range plotted. For
w
1.3
z a 0
I I
1.2
-
-~
UI
fully developed boundary layer flow one can cxpect that Sh is proportional to the ratio L over concentration boundary layer width 6~.Since ~ (Or)i/4 [9,19], no maximum should be observed. For comparison we have also plotted the corresponding KU curves (light lines) in fig. 14.
3,3.2. Inclination Now we will discuss the results from variations
1,.
1.0 -
0’
20’
40’
60’ 80’ 100’ 120’ INCLINATION ~s
140’ 160’ 180’
Fig. 15. Dependence of Sherwood number on inclination at AR 3. Note that maximum occurs between 9O0~I8(I0. i.e. not for the vertical, bottom-heated enclosure.
B. L. Markham, F. Rosenberger I
-
I
F
/
field on orientation change. In horizontal en-
~TT~T~
Pr =0.7 Sch =1.0 AR =10
more or less parallel to the side walls. On inclination, however, the convection roll becomes significantly tilted with respect to the side wall, reflecting the complex interaction of horizontal and vertical temperature (density) gradients in these systems. closures (4) = 90°) the core flow (fig. 6) proceeds we found (after considerable computational diffiFor the bottom-heated configuration (4) = 180°)
. -
6
-
8 ~ O~ 4 w I UI
culties) a complete change flow character. The asymmetric unicellular roll, in schematically depicted
3
-
0..r 0’
251
Diffusive— convective vapor transport
I 20’
I
40’
I I I J~,..I I 60’ 80’ 100’ 120’ 140’ i60’lBO’ INCLINATION
~
Fig. 16. Dependence of Sherwood number on inclination at AR = 10.
mere convective behavior in rectangular enclosures. Earlier numerical results show a maximum in the Nusselt number (thermal Sherwood number) at similar values [9,21,22]. Schinkel [9] has found for AR ~ 1, that the maximum Nusselt number is reached at 1200. For systems with adiabatic walls but AR’s and Or’s similar to ours [21,22] the Nusselt number has been found to reach its largest value at 130°.It appears that only Wirtz and Tseng [21] have considered a systern with conducting walls and an AR similar to ours. Their solutions indicate that the highest Nusselt number is reached at 1350. Suprisingly how—
—
ever, considering our work and that of Schinkel, these authors found that the Nusselt number stayed at anemometry constant studies value showed for larger also 4)’s. that Laser forVT the above Or maximum occur 900 <4) atthe < AR around transport 1800 the was maximum 4)vector = rates also 130° seen atwhich convection [14]. some in actual This angle trend velocities between experifor ments based on (ref. an [7], effective fig. 2). Or Hence, in transport simply aDoppler analysis projectiona and of gravity onto the main transport direction is used (eq. (55) in ref. [7]) cannot be expected to yield realistic results. The complexity of the Sh results are well reflected in the changes that occur in the velocity
in fig. 1, gave way to a double roll with symmetry about the middle plane y = h/2: upward flow in the middle of the enclosure and downward recirculation along the side walls. This is illustrated in fig. 17 by a U-velocity profile at x = L/2 for Or = 1.78 X 10”. Note that the area of the positive velocity exceeds those of the negative (downward)
regions. This is due to the net transport of A across (up) the cavity. For 4) = 00, i.e. for top-heating, we found also symmetry of the velocity profile with respect to the middle plane y = h/2, as depicted in fig. 17. This curve clearly shows the presence of convection even for top heating. Diffusive—advective mass transport from the hotter source tends to raise the temperature along the midplane more than closer to the (thermally stabilizing) side walls. Hence, horizontal temperature gradients and expansive convection develops, which leads to the character-
~
‘~
~
0,J _____ 0.6 0.4 0.2 0 -0.1
11” 0
I 0.1
~
180 Or Pr—0.7 =1.78x10’ Sch ~R =1.0 =10 0.2 0.3
DIMENSIONLESS VELOCITY
Fig. 17. Profile of velocity in transport direction in plane x = L/2 at 4 = 0 and 1800, showing symmetry about plane y = h/2.
252
B. L. Markham. F. Ro.senherger
/
istic dip in the velocity profile in fig. 17 (4) = 0). This we have discussed before in more detail in connection with solutal convection [10,23]. Of course, in treatments with uncoupled mass, momentum and energy transport, such insight cannot be obtained. Here we should also point out that a velocity profile similar to the 4) = 0 case of fig. 17 has been observed previously [18]. However, the authors interpreted this profile merely in terms of a diffusive deformation of a parabolic flow profile rather than the convective flow that results from the horizontal density gradients in such a configuration. The observed flow symmetry with respect to = h/2 is of some fundamental interest. For mere convection in vertical cylinders, rolls at not too high Or are asymmetric. In our earlier work on VT in vertical cylinders, however, we had assumed axisymmetric flows (for numerical simplicity’s sake) and argued that these will prevail due to the (symmetric) net transport in the system [101. The above findings give limited (i.e. 2D) support to this assumption. Finally, we have investigated the growth flux non-uniformities that arise in inclined cavities. In general, for 4) < 90° the non-uniformity is somewhat reduced as compared to fig. 7. For 4)> 90° it is first slightly increased up to the maximum transport angle, 4)max and then reduced. In fig. 18 we compare the growth flux distribution for AR = 10, Or = 1.78 X iO~and 4)fu~ = 130° (see fig. 16) with
~
1.0
\
Or =i.78xi0’ 0.8
/
AR =10 06
~/ ,,/
-
90~~>/
/~
0.4
0.2
/ 7 /
-
f 0
\ )
~ ,/
~
_•_7~
0.6
0.8
1.0
lower (upper) half of the enclosure. The remaining part changes the flow direction by dropping (rising) at the side wall over a considerably distance from the respective end. Further evidence of the fundamental shortcon3ings of 2D discussions comes from current 3D
below our model predictions.
1.21.4
DIMENSIONLESS GROWTH FLUX
Fig. 18 Dimensionless growth flux distribution for and 130°.
In order to compare our predictions with cxperimental results we have performed a detailed literature search, There are, however, only few VT experiments that have been characterized to an extent that encourages a comparison [6,7]. Even with these experiments we can only obtain semiquantitative correlation. One could list numerous parameters that are not accounted for even in our computationally complex model and that may he important in reality. Among these are the neglect of radiative heat transfer that tends to reduce the effective Or, the often more complicated interface geometries obtained in reality and the neglect of the Soret effect. Yet, we believe that the lack of quantitative agreement has more fundamental rcasons, Recent laser Doppler anemometry studies of the convective velocity field in cylinders performed in our laboratory have revealed rather unexpected 3D flow features [13]. It was found that only part of the vapor flows along the cold (hot) thermode. i.e. the “crystal” and “source” in fig. I, into the
functions of figs. 3 and 4 ought to be shifted to the right. This trend is supported by the observation that most experimental results for Sh appear to lie
,/__
0.4
3.4. Comparison with experimenl.v
of the Fluid Dynamics Oroup at the University of Marseille and our group [15,16]. These efforts show that, in general, 2D systems enter the various fluid dynamics regimes. (such as boundary layer flow) at lower Or values than the corresponding 3D systems. This implies, for instance, that the Sh(Gr) .
/
02
the distribution obtained for 4) = 90°. Besides the increase in transport one notes also a shift in the location of the flux maximum as a result of the before mentioned tilt of the convection roll.
simulation that is being performed in a joint effort
\ \
-
Diffusive — conI’ectil’e vapor tra/isporl
~
=
900
Why then strive for improved 2D models, as we have attempted in this work? Full-fledged 3D treatments of the three coupled transport equa-
B.L. Markham. F. Rosenberger
/
tions for mass, momentum and energy (including
concentration-induced density gradients) for realistic boundary conditions will remain beyond reach for most workers for some time to come. Hence, 2D models are currently the only theoretical means of obtaining insight on practically important trends that can lead to improvements of the transport conditions.
4. Summary Our numerical solutions for horizontal VT agree well with predictions of the KU model at low Grashof numbers and large aspect ratios. At medium and high tend to level off. Or, Our KU’s modelSherwood yields a numbers realistic increase which asymptotically approaches a (Gr)°25 dependence, in agreement with boundary layer theory. At low aspect ratios and moderate products of Sch and Or, KU’s analysis tends to overestimate the Sh. In both approaches a maxi-
Diffusive— convective vapor transport
NSO-1534 and by the National Science Foundation under grant DMR-79-13183 are gratefully
acknowledged. Stimulating discussions with P. Bontoux are also very much appreciated.
Nomenclature A AR
Gaseous growth component Aspect ratio (L/h)
B C 1,
Gaseous inert component
DAB g Or h k L p 1 P Pc
mum in Sh is reached at smaller aspect ratios, the exact value of AR slowly growing with Or as long as full boundary layer behavior is not achieved. The Schmidt and Prandtl numbers of the vapor were shown to have significant influence on transport, effects that cannot realistically be deduced from the simplifying KU model. We have also investigated the effect of orientation of the enclosure on Sh. Transport becomes maximum for an angle of about 300_400 between the main transport direction and the horizontal with the hotter end below. The exact value of this angle, however, must be expected to depend on Or and Ar. Beyond the overall transport predictions, our analysis yields also solutions for temperature and growth flux distributions. These show practically important non-uniformities, and shed light on the physical origin of shortcomings of earlier transport models.
Specific heat at constant pressure per unit mass Binary diffusion coefficient Acceleration due to gravity 3~T/~2) Orashof Height ofnumber cavity ($gh Thermal conductivity Length of cavity Partial pressure of ith (A or B) component Total pressure Péclet number (ln[ ~vB( L )/p ~~(0)], see ref.
[3]) Pr Sch Sh
T L~T T
Prandtl number (v/K) Schmidt number (v/DAB) Sherwood number (ratio of total interfacial transport fluxes in 2D and corresponding 1D model) Temperature T~— T~
Normalized temperature [(T T )/( T~ 1’)] —
7~
T~ u
U v V
W~ x X y
Acknowledgments Y Support of this work by the National Aeronautics and Space Administration under grant
253
—
Temperature at crystal Temperature at source Component of velocity in x-direction Dimensionless x-component of velocity (uh/v) Component of velocity in y-direction Dimensionless y-component of velocity (vh/v) Weight fraction of I th (A or B) component Coordinate in main transport direction Dimensionless x-coordinate (x/L) Coordinate perpendicular to main transport direction Dimensionless y-coordinate (y/h)
254
1k!.. Markham. F Rosenherger
/
DiJfuso’e-— ionIeelIvr’ vapor trailsport
Greek letters
Enclosures (Dutch
Efficiency Bureau-- Pijnacker.
t)elft.
1980).
,8
Thermal expansion coefficient (ideal gas:$ = 1 /T / . . . Thermal diffusivity Viscosity Kinematic viscosity Mass density
t~
p. v
p
12]
Crystal Growth 51(1981) 426. 111 R. Jackson. Transport iii Porous Amsterdam, 1977) eli. 4.
Catalysts
Elsesier,
112) B.S. Jhaveri. B.L. Markham and F. Rosenherger. (heni. Eng. Comntun. 13(1981) 65.
[131(LII.
Schirok~and F. Rosenherger. Intern. .1. Heat Mass Fransfer 27 (1984) 587, [14] F. Rosenberger and AC. 1 Iurford. in preparation. 115] C. Smutek. P. Bontoux, B. Rou\. (iii. Schiroky, AC. Hurford. F. Rosenherger and U. de VahI Davis, Nunieri-
References )l] H. Wiedemeier. F.C. Klaessig. E.A, Irene and Si. Wev. Crystal Growth 31 (1975) 36.
101 BE. Markham, [).W. Greenwell and F. Rosenhcrgcr, J.
j.
H. Wiedemeier. J. Sadeek, F.C. Klaessig, M. Norek and R. Santandrea, J. Electrochem. Soc. 124 (1977) 1095.
cal Heat Transfer, submitted.
116] P. Bontoux. B. Roux, G.H. Schiroky, B.L. Markham and F. Rosenberger, Intern. J. 1-leat Mass Transfer. submitted. )17[
A.P. Gosman and W.M. Pun. Report No. UTS/74/2. Impenal College of Science and Technology, London,
[3] D.W. Greenwell, B.L. Markham and F. Rosenherger. J. Crystal Growth 51(1981) 413. [4] A. Solan and S. Ostraeh, Convection Effects in Crystal Growth by Closed-Tube Chemical Vapor Transport, in: Preparation and Properties of Solid State Materials. Vol.
1973. [181 G. Wyatt. W.E. Ihele and ER. Eckert. Intern. .1. Heat Mass Transfer 14(1973) 513. [19) A.E. Gill, J. Fluid Mech. 26 (1966) 515. [20] F. Rosenberger, Fundamentals of Crystal Growth I
4, Ed. W.R. Wilcox (Dekker, New York. 1979) pp. 63--I 10. [5] K. Klosse and P. Vllersnia, J. Crystal Growth 18 (1973) 167. [6) C’. Paorici and C. Pelosi. J. (‘rystal Growth 35(1976) 65.
(Springer, New York. 1979) p. 300. 1211 R.A. Wirtz and W.F. Tseng, Natural Convection Across Tilted Rectangular Enclosures of Small Aspect Ratio, in: Natural Convection in Fnclosures—EITD. Vol. 8. Eds. J.
[71 D.
Chandra and H. Wiedemeier, J. Crystal Growth 57 (1982) 159. [8] B.S. Jhaveri arid F. Rosenberger. J. Crystal Growth 57
Catton and K.E. Torrance (ASME, New York. 198W pp. 47—54. 22] iC. Grondin, Phi) Thesis. Université de Prosence.
(1982) 57. [9] W. Schinkel, Natural Convection in Inclined Air-Filled
Marseille (1978). [23) B.L. Markham and F. Rosenberger, Chem. Eng. (‘ornniun. 5(1980) 287.