# Numerical simulations of buoyancy driven flows in cylinders and cavities for vapour crystal growth

Patrick Bontoux, Frédéric Elie, Claude Smutek, Guy Pierre Extremet, Anthony Randriamampianina, Emilia Crespo Del Arco, Hubert Branger, Bernard Roux

► To cite this version:

Patrick Bontoux, Frédéric Elie, Claude Smutek, Guy Pierre Extremet, Anthony Randriamampianina, et al.. Numerical simulations of buoyancy driven flows in cylinders and cavities for vapour crystal growth. Manuel G. Velarde. Physicochemical hydrodynamics. Interfacial Phenomena, 174, Kluwer Academic Publisher, pp.963-995, 1989, NATO ASI Series, 978-0-306-42905-7. 10.1007/978-1-4613-0707-5_68. hal-00229875

HAL Id: hal-00229875

https://hal.science/hal-00229875v1

Submitted on 19 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
NUMERICAL SIMULATIONS OF BUOY ANCY DRIVEN FLOWS IN CYLINDERS
AND CAVITIES FOR VAPOUR CRYSTAL GROWTH

P. Bontoux, F. Elie, C. Smutek, G. P. Extremet',
A. Randriamampianina, E. Crespo *, H. Branger, and B. Roux

I.M.F.M., rue Honnorat, Marseille, France

INTRODUCTION

The paper reports about numerical works carried out by the
authors at I.M.F.M. in the group of Numerical Fluid Mechanics. The
common objectives of these works are related to the analysis of fluid
dynamics in closed tube methods employed in crystal growth from
vapours. Several numerical methods are considered and the solutions
are discussed with respect to the results obtained with different
approaches of the problem as the asymptotic theories, the stability
analyses and the experiments. Some groups were associated in relation
with these approaches and also the development of the numerical
methods, finite difference and spectral methods. They are the groups of
Prof. G. De Vahl Davis at University of New South Wales-Australia, Prof.
F. Rosenberger at University of Utah-U.S.A. (with Dr. G.H. Schiroky, B.L.
Markham, A.C. Hurford..), Prof. R. Sani at University of Colorado-
Boulder-U.S.A. (with G. Hardin), Prof. M.G. Velarde at U.N.E.D.-Madrid-
Spain (with E. Crespo), Dr. R. Peyret at Université de Nice-France (with
Dr. J. Ouazzani and J.M. Lacroix).

The closed ampoule vapour transport techniques are used for the
preparation of single crystal and very useful for research in crystal
growth (however costly for production where open tube methods can be
generally preferred). The ampoule is usually cylindrical and heated
differentially inside a furnace which produces the thermal conditions for
vaporization at the source and condensation at the sink (see Omaly et al.,
1981). The process involves very complex physical mechanisms, i.e. the
sublimation of a material at the hot source, the transport of the volatile
component A across the ampoule through other components B and the
crystal growth at the cold end of the tube (see Rosenberger and co-
workers, 1979 to 1984). The convection is the principal fluid mechanism
and the transport can also be enhanced by tilting the tube with respect
to the gravity field (Pamplin, 1980). Depending on the values of the

<footnote>*permanent address: Dpto. Fis. Fund., U.N.E.D., Madrid, Spain</footnote>
physical parameters (Rayleigh numbers...), the aspect ratio and the inclination angle, many complex flows can develop in the ampoule which can be highly or moderately three-dimensional or simply two-dimensional as in axisymmetric situations. The understanding of convection in the nutrient by predictions of velocity scales and flow patterns is very important for the researchers because it affects the growth behaviour of the crystal (Carruthers, 1977, Olson and Rosenberger, 1979, Ostrach and co-workers, 1979, 1981, 1982).

The mathematical model is given by the Navier-Stokes equations with the Boussinesq approximation. The three-dimensional simulation of thermal convection is made with a finite-difference method. Both hermitian finite differences and spectral approximations were considered for the two-dimensional simulation of physical vapour transport with a simple analytical law or the Fick's law to modelize the mass flux at the interfaces. The characteristics of the numerical methods are shortly given. The typical flows obtained in horizontal, vertical and inclined cylinders are presented. An assessment of the validity of the two-dimensional approximation is discussed with respect to the experimental, analytical and numerical results in the symmetry plane of a horizontal cylinder. The fluid dynamics of crystal growth from vapours is studied when the thermal and solutal Rayleigh numbers vary.

## PHYSICAL MODELS AND GEOMETRIES

Crystal growth by physical vapour transport (PVT) in a closed tube is experimentally simple. It consists in the vaporization of a material A at a source and its condensation to a crystal generally held at a "colder" temperature than the source (Greenwell et al., 1981, Markham et al., 1981). The transport of the gaseous component A through (inert) components B across the ampoule is governed by the competition between the advection-diffusion driven by the interfacial mass flux and the convection generated by the density gradients in the changes of temperature and concentration.

The cylindrical geometry is depicted in Fig.1 with the associate rectangular cavity model. The reference frames are given for both geometries with the velocity components $u, v, w$ defined in the radial ($r$ or $x$), azimuthal ($\phi$) and axial ($z$) directions, respectively. The axial length is $L$, $R$ is the radius and $D$ is the height ($D=2R$). The aspect ratio is $a=L/2R=L/D$. The enclosure is filled with a binary mixture of gaseous components A and B (molecular weights $M_A$ and $M_B$). At the source and the crystal, the solid S consists of A only. The interfaces (of length D) are flat, stationary and held at constant different temperatures (hot at the source $T_2$ ($z=L$), cold at the sink $T_1$ ($z=0$)) and concentrations ($S_{A2} \ge S_{A1}$). The growth rate is assumed to be limited by the mass transfer of A in the vapour phase. The flux of the (inert) component B vanishes at the interfaces (active walls). The passive side walls are impermeable to both species.
<figure><img src="image_2.png" /><figcaption>Fig.1. Differentially heated cylinder and Fig.2. Temperature profiles along the cavity at inclination with the gravity (axial) walls: 1-G model (a), 2-Z model vector: definition sketch for dimensions, (b) and 3-Z model (c). position coordinates and velocity components.</figcaption></figure>

<figure><img src="image_3.png" /></figure>

Depending on the applications, various kinds of temperature profiles can be produced at the axial walls by the furnace in order to enhance the growth process (Omaly et al., 1981, 1983, Launay, 1982, Brisson, 1981, Extrémet et al., 1986). They are presented in Fig.2. The simplest profile is linear and corresponds to the one-gradient model (abbreviated 1-G). In the general case one or more gradient can be inserted between isothermal zones (multizone models, with two and three zones abbreviated 2-Z and 3-Z). The conditions at the side walls are expressed as:

$$T = T_1 + \Delta T \Theta(z) \quad (1)$$

where $\Theta(z)=z/L$ with $\Delta T=T_2-T_1$, when $0 \le z \le L$ for the 1-G model (Fig.2a); $\Theta(z)=0$ when $0 \le z \le z_1$ and $\Theta(z)=1$ when $z_2 \le z \le L$ for the 2-Z model (Fig.2b); $\Theta(z)=0$ when $0 \le z \le z_1$, $\Theta(z)=\Theta_m=(T_m-T_1)/\Delta T$ when $z_2 \le z \le z_3$, $\Theta(z)=1$ when $z_4 \le z \le L$ for the 3-Z model (Fig.2c); in the 2-Z and 3-Z models the constant gradient zones (of widths $\delta_i=(L_{\Delta T})_i/L$) are inserted between the isothermal zones and $T_1 \le T_m \le T_2$.

The axis of the ampoule will not be always perpendicular to the gravity field, that is in microgravity environment but also on earth in the applications when the inclination is used to enhance convection (Pamplin, 1980). As defined in Fig.1, the inclination angle, $\Omega$, is referred to the vertical. The inclination $\Omega=0^\circ$ will correspond to the vertical cavity heated from below.
MATHEMATICAL MODELS

The governing system is given by the Navier Stokes, energy and species equations with the Boussinesq approximation and the linear Fick's law.

The Boussinesq Approximation and the Linear Fick's Law

The density is linearly related to the variation of temperature, T, and mass fraction, $S_A = \rho_A/\rho$, of component A ($S_A+S_B=1$) with the Boussinesq approximation (Joseph, 1976) as follows :

$$ \rho = \rho_0 \{ 1 - \beta (T - T_0) - \alpha_A (S_A - S_{A0}) \} \quad (2) $$

where $\beta$ is the thermal expansion factor ($1/T_0$ for ideal gas), and $\alpha_A = M_0(M_B - M_A)/M_A M_B$ is the solutal expansion factor expressed in terms of the molecular weights (subscript 0 refers to average quantities).

The diffusive flux of the components with respect to the mass average velocity, $V=S_A V_A + S_B V_B$, is given by the linear Fick's law as (see Rosenberger, 1979):

$$ S_A(V_A - V) = -D_{AB} \nabla S_A \quad (3) $$

where $D_{AB}$ is the binary diffusivity.

The Two-Dimensional Model

The governing equations are written using the dimensionless vorticity, $\zeta$, streamfunction, $\psi$, temperature, $\theta=2(T-T_0)/(T_2-T_1)$, and mass fraction, $S=2(S_A-S_{A0})/(S_{A2}-S_{A1})$, as :

$$ \frac{\partial \zeta}{\partial t} + \eta_i V \cdot \nabla \zeta = \eta_d \Delta \zeta + \eta_{tb} g_0 \cdot \nabla \theta + \eta_{sb} g_0 \cdot \nabla S \quad (4) $$

$$ \frac{\partial \theta}{\partial t} + \lambda_i V \cdot \nabla \theta = \lambda_d \Delta \theta \quad (5) $$

$$ \frac{\partial S}{\partial t} + \sigma_i V \cdot \nabla S = \sigma_d \Delta S \quad (6) $$

$$ \Delta \psi = \zeta \quad (7) $$

where $\nabla$ and $\Delta$ are the gradient and Laplacian operators, $V=(u=\partial\psi/\partial z, w=-\partial\psi/\partial x)$ and $g_0$ is the unity gravity vector. Depending on the choice of the scaling factors for the velocities and space variables, the coefficients $\eta, \lambda, \sigma$ are expressed in terms of the physical parameters (De Vahl Davis, 1986): the Prandtl number, $Pr=\nu/\kappa$, the Lewis number, $Le=\kappa/D_{AB}$, or the Schmidt number, $Sc=\nu/D_{AB}$, the thermal and solutal Rayleigh numbers, $Ra_T=\beta g \Delta T (H)^3/\nu \kappa$ and $Ra_S=\alpha_A g \Delta S_A (H)^3/\nu D_{AB}$, or
the Grashof numbers, Gr=Ra/Pr, where $\Delta S_A=S_{A2}-S_{A1}$ and H=D/2 (or R in what will concern the cylindrical cavity). With $\kappa/H$ ($\kappa$ is the thermal diffusivity) as scaling factor for the velocity, the dimensionless coefficients in system (4-7) are $\eta_i=\lambda_i=\sigma_i=1$, $\eta_d=\text{Pr}$, $\lambda_d=1$, $\sigma_d=1/\text{Le}$, $\eta_{tb}=\text{Ra}_T\text{Pr}$ and $\eta_{sb}=\text{Ra}_S\text{Pr}/\text{Le}$. The common scaling for the space variables in the following results will be based on the mid-height D/2 of the cavity (respectively R for the cylinder).

The mass fraction at the side walls is governed by the no-flux condition for both species, $\partial S/\partial x=0$. For the velocity the usual no-slip and no-permeability conditions are used (Jhaveri et al., 1981, 1982).

Two models are considered to analyse the effect of the interfacial mass flux of component A at the active walls:

- one is derived directly from the Fick's law (3) (see Rosenberger and co-workers, 1981 to 1984); the w-profile is then:

$$w_{\text{interface}} = (\partial S / \partial z) \xi(x) / \text{Le} (E - S) \quad (8)$$

where E is a dimensionless number characterizing the mass fractions at the source and the sink (and then related to the molecular weights $M_A$ and $M_B$ and partial pressure conditions at sublimation and deposition, see Elie, 1984):

$$E = 2 (1 - S_{A0}) / \Delta S_A . \quad (9)$$

The transport of components A and B is computed from the mass average solution using a second number F, defined as $F+E=2/\Delta S_A$. The function $\xi(x)$ is artificially introduced in order to regularize the solution near the corners. Its constant is also adjusted at the sink in order to ensure the mass conservation of species A (see Ouazzani, 1984, Elie, 1984).

- an analytical model is also considered as a perturbation of the convective flow and following the quartic w-velocity profile (Extrémet et al., 1986)

$$w_{\text{interface}} = v_d q [x(x-1)]^2 \quad (10)$$

where q is related to the maximal velocity generated by buoyancy in the core (see following relation (14)) and $v_d$ is a proportionality factor suitably chosen to correspond to realistic conditions.

## The Three-Dimensional Model

The governing system is considered with the velocity (V)-vorticity ($\zeta=\nabla \times V$) formulation where $V(u,v,w)$ (see definitions in Fig.1). With R and $\kappa/R$ as scaling factors for the coordinates and the velocities and
with θ as dimensionless temperature, the governing equations for
thermal convection are

$$
\frac{\partial \zeta}{\partial t} = \nabla \times (V \times \zeta) - 0.5 \mathrm{Ra}_T \mathrm{Pr} \nabla \times (\theta g_0) - \mathrm{Pr} \nabla \times (\nabla \times \zeta) \quad (11)
$$

$$
\frac{\partial \theta}{\partial t} = - \nabla \cdot (V \cdot \theta) + \nabla^2 \theta \quad (12)
$$

$$
\nabla^2 V = - \nabla \times \zeta + \nabla (\nabla \cdot V) \approx - \nabla \times \zeta \quad (13)
$$

The thermal boundary conditions correspond to the 1-G model. In addition to the usual no-slip conditions, the conditions on the vorticity at the rigid walls are derived as : $\zeta_r = -\partial v/\partial z$, $\zeta_\phi = \partial u/\partial z$, $\zeta_z = 0$ at the end walls and $\zeta_r = 0$, $\zeta_\phi = -\partial w/\partial r$, $\zeta_z = \partial v/\partial r$ at the side walls.

ANALYTICAL SOLUTIONS

The basic flow in long horizontal cavities or cylinders filled with a
mono-component fluid corresponds to two counter-flows coming from
the hot zone in the upper part, and from the cold zone in the lower part.
Analytical approximations for the axial velocity profile are available for
the core of rectangular cavities and cylinders when they are
differentially heated using a 1-G temperature model (see Bontoux et al.,
1986). These approximations are based on temperature gradient
parameters derived from experiments and analyses. At low Ra<sub>T</sub>, the
horizontal temperature gradient is constant and generates the main
buoyancy forces, which corresponds to the core driven regime, denoted
CDR (see Cormack et al., 1974, Imberger, 1974). At larger Ra<sub>T</sub>, the main
buoyancy forces are located in the end regions (boundary-layer driven
regime, denoted BLDR). In these cases the functional laws available from
Cormack et al. (1974), Bejan and Tien (1978) and Hart (1983) can be
used to predict the flow profiles. The analysis was also recently
extended for cavities in the case of multizone models by Extrémet et al.
(1986).

1-G Thermal Model for Cavities and Cylinders

Parallel flow solutions were proposed by Klosse and Ullersma
(1973), Birikh (1966), Hart (1972), and Bejan and Tien (1978). Elaborate
approximations, including some interaction with the end regions and
allowing for secondary flows in cavities, were given by Cormack et al.
(1974), Bejan and Tien (1978) and Shih (1981) (see also Ostrach et al.,
1980, 1982). These core solutions are still limited to the S-shaped
profile for the horizontal velocity component. Further improvements,
which remove this limit for higher Ra<sub>T</sub> conditions, are due to Tichy and
Gadgil (1982). For cylinders, Schiroky and Rosenberger (1984) proposed
a third-order power series in Ra<sub>T</sub> which well predicts the shift of the
maximum in the core velocity profile towards the wall for the low Ra<sub>T</sub>-
range of the BLDR.
The analytical expressions for the core flow (horizontal velocity) depend on the axial dimensionless temperature gradient $k_1=0.5(\partial\theta/\partial z)$. They are expressed as follows in terms of $k_1 Ra_T$:

- for cavities:

$$w^{2D} = (k_1 Ra_T) (x^2-1)x/6 \quad (14)$$

- for cylinders:

$$w^{3D} = (k_1 Ra_T) \cos\phi (r^2-1)r/8 \quad (15a)$$

$$u^{3D} = -(k_1 Ra_T)^2 \cos 2\phi (2r^6 - 15r^4 + 24r^2 - 11)r/184320 \quad (15b)$$

$$v^{3D} = -(k_1 Ra_T)^2 \sin 2\phi (8r^6 - 45r^4 + 48r^2 - 11)r/184320 \quad (15c)$$

In the CDR the value of $k_1$ is $(2a)^{-1}$ and the flow is mainly axial. For the CDR up to the beginning of the BLDR, Cormack et al. (1974) derived an expression for $k_1$ in terms of $Ra_T^2$ and a. For the fully developed BLDR, Bejan and Tien (1978) proposed a $Ra_T^{-3/5}$ dependence. More recently, Hart (1983) derived from Cormack et al's works an expression for $k_1$ that is valid for the whole CDR and BLDR range. The expression is implicit but, for large values of $Ra_T$, one recovers the explicit relation

$$k_1 = (1024Q)^{-1/3} Ra_T^{-2/3} \quad (16)$$

which is independent of the aspect ratio a (Q is a constant parameter), as Bejan and Tien's relation. For cylinders Schiroky (1982) proposed an extension of the end-integral method used by Bejan and Tien (1978) for the prediction of $k_1$ (see also in Bontoux et al., 1986). When $Ra_T$ is increased, three-dimensional flow structures superimpose to the basic flow as described by relations (15b-c).

## Multizone Thermal Models for Cavities

With a multizone temperature profile, the analysis was developed by Extrémet et al. (1986) for the conduction regime (CDR). The temperature gradient in the middle of the i'th temperature-gradient zone is not only connected to the aspect ratio a as in the 1-G model but also depends on the reduced width of the temperature gradient $\sigma_i$. The functional laws for $k_1$ are empirically determined from computations (see Fig.3) when the wall gradient is inserted between two isothermal regions (2-Z model) as

$$2k_1=1 \quad \text{at} \quad a\sigma_i<1 \quad (17a)$$

$$2a\sigma_i k_1 = 1 \quad \text{at} \quad a\sigma_i > 1 \quad (17b)$$
<figure><img src="image_4.png" /><figcaption>Fig.3. Variation of the temperature gradient k₁ in 2-Z models with 10≥a≥1 and 1≥γ≥0.05 for conduction regime. From Extrémet et al., 1986.</figcaption></figure>

The maxima of the velocity profile (14) can be expressed by:

$$|w^{2D}_{max}| = A/\sqrt{3} \qquad (18)$$

with A= Ra_T k_1 β_i δ_i and where β_i is a discrepancy factor determined from comparisons with computed solutions (Extrémet et al., 1986). β_i was found to vary between 1 and 0.4 for a wide range of aspect ratios, a, and reduced widths of the temperature gradient, δ_i. The factor δ_i is 1 for the 2-Z model and is proportional to the local reduced temperature step, as Θ_m=(T_m-T_1)/ΔT and (1-Θ_m) for the 3-Z model (see Extrémet et al., 1986).

## Binary Mixture

The parallel flow solution is extended to the case of thermal and solutal convection, as follows, where kS₁ is related to the axial concentration gradient, kS₁=0.5(∂S/∂z):

$$w^{2D} = (k_1 Ra_T \beta_i \delta_i + k^{S_1} Ra_S / Le) (x^2 - 1)x/6 \qquad (19)$$

The motion corresponds again to two counterflows, but which direction depends on the sign of (k₁Ra_Tβ_iδ_i + kS₁Ra_S/Le). For the 1-G model and at low Ra_T, the horizontal gradients of temperature and concentration are constant in the core. Both these gradients are then nearly (2a)⁻¹. At the large Ra_T, the functional laws available from Bejan and Tien (1978) and Hart (1983) can be used as a first approximation for a binary fluid. Also, the solution (19) was discussed with respect to computed solutions by Extrémet et al. (1986) for multizone models.
The solution of governing systems (4-7) and (11-13) was made by using different numerical methods. Classical finite differences known for their "robustness" were used for the computation of the complex three-dimensional flow patterns in cylinders. In cavities filled with a binary mixture higher-order accurate methods were considered as the hermitian finite-difference method and the spectral method based on Chebyshev-polynomial expansion of the variables. Their characteristics are briefly recalled with the details corresponding to the present applications.

Classical Finite-Difference Methods

The approximation is local and based on the classical explicit relations between the derivatives and the variables over three discretizing points. With a constant step size in space, Δz, the relations are second-order accurate and written as:

$$ (\partial f / \partial z)_k = (f_{k+1} - f_{k-1}) / 2 \Delta z \quad (20a) $$

$$ (\partial^2 f / \partial z^2)_k = (f_{k+1} - 2f_k + f_{k-1}) / \Delta z^2 \quad (20b) $$

The details of the method developed for the three-dimensional convection problem are reported in Leong and De Vahl Davis (1979), Leong (1983) and Smutek et al. (1983, 1985). The method is based on centered finite differences with an uniform mesh composed of L x M x N discretizing points in r-, φ- and z-directions as $r_i = (i - 1/2) \Delta r$ with i=1 to L, $\phi_j = (j - 1) \Delta \phi$ with j=1 to M and $z_k = (k - 1) \Delta z$ with k=1 to N, where $\Delta r = 1/(L - 1/2)$, $\Delta \phi = 2\pi/M$, $\Delta z = 2a/(N - 1)$. No mesh points are located on the axis to avoid singularity problems (De Vahl Davis, 1979).

The continuity equation, $\nabla \cdot V = 0$, is not automatically satisfied. Then, it is used to check the accuracy resulting from a given mesh size (see Bontoux et al., 1986) The solution of Poisson velocity equations (13) is made with the Fourier series direct method (Le Bail, 1972) and using the FFT algorithm (Cooley and Tukey, 1965). The advancement in time is based on Samarskii-Andreyev ADI scheme (Mallinson and De Vahl Davis, 1973, Samuels and Churchill, 1967, Peaceman and Rachford, 1955).

The convergence of the solution towards a steady state with a 9x32x33 mesh can be obtained with less than 100 iterations up to $Ra_T \approx 4000$ and starting from initial conditions derived from relation (15) at $a=5$ and $Pr=0.73$. After the transition to the boundary layer driven regime (BLDR), convergence is generally obtained by incrementing $Ra_T$. For the most severe conditions at $Ra_T=18,720$ and when the axis is tilted from the horizontal to nearly the vertical ($90^\circ \ge \Omega \ge 20^\circ$) the convergence
is achieved after about 500 to 3500 iterations using a 9x32x65-mesh and varying progressively Ω by steps of 10° to 15°. The computation cost per iteration is about 0.32, 0.60 and 1.19 sec on CRAY1/S with 9x32x33, 65 and 129-meshes, respectively.

## Hermitian Finite-Difference Methods

The method considers implicit fourth-order accurate relations between the derivatives and the variable at three discretizing points (Krause and co-workers, 1972-76, Hirsh, 1975). The derivatives are considered as additional unknowns of the problem but can be eliminated in the case of linear solvers as in Poisson's equation (Adam, 1977). The implicit relations are:

$$ (\partial f/\partial z)_{k+1} + 4(\partial f/\partial z)_k + (\partial f/\partial z)_{k-1} = 3(f_{k+1} - f_{k-1})/\Delta z \quad (21a) $$

$$ (\partial^2 f/\partial z^2)_{k+1} + 10(\partial^2 f/\partial z^2)_k + (\partial^2 f/\partial z^2)_{k-1} = 12(f_{k+1} - 2f_k + f_{k-1})/\Delta z^2 \quad (21b) $$

The method was applied to the two-dimensional governing system, using for the closure of the system additional relations at the boundaries as the Padé approximant (Roux et al., 1978, Bontoux et al., 1978, Peyret and Taylor, 1982). With the vorticity and streamfunction formulation a third-order accurate relation is considered for the vorticity at the walls. A mixed method was used here and based on classical finite differences for the transport equations combined with the hermitian method for Poisson's equation and Hirsh's relation for the wall vorticity (Loc and Daube, 1978, Roux et al., 1980). The integration is made with an A.D.I. false-transient method (see Peaceman and Rashford, 1955, Mallinson and De Vahl Davis, 1973). The false-transient method is used which introduces relaxation factor to improve the convergence process.

For a converged solution obtained after 240 iterations with 21x81 discretizing points the computations take 3.15 seconds on CRAY1/S computer.

## Spectral Methods

The spectral methods are based on the method of the weighted residuals and on trial-function expansions of the variables. Their use has become very popular as numerical methods after pioneering works by Orszag (1971) (see Gottlieb and Orszag, 1977) and the development of vector computers. The approximation gives highly accurate derivations in the spectral plane. The method is based on powerful vectorized algorithms as Fast Fourier Transforms (Temperton, 1983) and Poisson-Helmholtz solvers (Haidvogel and Zang, 1979, Bondet de la Bernardie, 1980).

In the case of Chebyshev-polynomial expansions the variable and its derivatives are expressed as:
$$f \approx \sum_{n=0}^{N} a_n T_n(x) \qquad (22a)$$

$$(\partial f / \partial z) \approx \sum_{n=0}^{N} a_n T_n(x) \qquad (22b)$$

$$ \begin{array}{l} \text{(1) N} \\ \text{with } a_n = (2/c_n) \sum_{\substack{p=n+1 \\ p+n \text{ odd}}}^{N} p a_p \end{array} $$ $$ (\partial^2 f / \partial z^2) \approx \sum_{n=0}^{N} a_n T_n(x) \qquad (22c) $$ $$ \text{(2) N} \\ \text{with } a_n = (1/c_n) \sum_{\substack{p=n+2 \\ p+n \text{ even}}}^{N} p^2 a_p $$

The solution of system (4-7) is made with the Tau-Chebyshev spectral method using Chebyshev polynomial expansions both in x- and z-directions (see Bondet de la Bernardie, 1980, Elie et al., 1983, Randriamampianina, 1984) as

$$(\zeta, \theta, \psi, S) \approx \sum_{n=0}^{N} \sum_{m=0}^{M} (a, b, c, d)_{nm} T_n(x) T_m(z) \qquad (23)$$

It results in a set of ordinary differential equations in the spectral space for the transport equations as:

$$ \begin{array}{l} \displaystyle \frac{da_{nm}}{dt} + \eta_i ( e_{nm} + e_{nm}^* ) - \eta_d ( a_{nm} + a_{nm}^* ) - \eta_{tb} b^{*nm} - \eta_{sb} d^{*nm} = 0 \qquad (24) \\ \displaystyle \qquad 0 \le n \le N-2 \quad 0 \le m \le M-2 \quad \text{and} \quad 0 \le \zeta \le \theta \le \psi \quad \text{and} \quad S \text{ (diffusion and source terms).} \end{array} $$

and for the Poisson equation as

$$ c_{nm} + c_{nm}^{*} = a_{nm}, \qquad (25) $$

with $0 \le n \le N-2$ and $0 \le m \le M-2$ and where $b^{*}_{nm}$ and $d^{*}_{nm}$ are expressed similarly as $b^{*}_{nm} \sim \cos\Omega b^{10}_{nm} + \sin\Omega b^{01}_{nm}$. The superscripts ij for the variables a, b, c and d refer to the components of the ith x- and jth z-derivatives of $\zeta, \theta, \psi$ and S (diffusion and source terms). The components $e_{ij}, h_{ij}$ and $g_{ij}$ correspond to the non-linear convective
terms. The closure of the system with the conventional boundary conditions was detailed elsewhere (Bontoux et al., 1981, Elie, 1984, Randriamampianina, 1984).

Among the various FFT algorithms available for the evaluation of the non-linear terms (Singleton-IMSL, CRAY-SCILIB library, Lhomme et al., 1982, Temperton, 1983), the algorithm of Temperton is chosen as it saves computation time on vector computers and works for any numbers of Fourier components such that $2^{P345}$. Poisson equation (25) is solved using the Matrix Diagonalization Technique (Haidvogel and Zang, 1979, Bondet de la Bernardie, 1980).

Various schemes can be used for the integration in time (Adams-Bashforth, Adams-Bashforth-Crank-Nicolson schemes, Runge-Kutta, LSODA methods). The multistep LSODA method (Livermore Solver for Ordinary Differential Equations with Automatic method switching between non-stiff and stiff problems, ODEPACK) utilizes two predictor-corrector schemes adapted for non-stiff and stiff problems, respectively the Adams-Moulton scheme (AM) and the Backward-Differentiation Formula (BDF) (for details see Hindmarsh, 1976, 1982, Petzold, 1980).

With the LSODA method, no crucial stability problem was encountered as both the time step and order of integration are adapted automatically. The cost of LSODA is approximatively the same as the explicit AB method at small (NxM). When (NxM)>400 the iteration cost (as the memory required) increases strongly due to the internal generation of the jacobian which varies as 2(NxM) (see Randriamampianina et al., 1985).

## RESULTS

The results bring some insights about three points of interest concerning the prediction of vapour transport in closed tubes.

- The first point is related to the understanding of the complex three-dimensional flows which develop in cylinders when the buoyancy force is governed by the temperature differences and when the inclination varies. The finite difference solutions are analysed with respect to experimental and theoretical results when they exist.

- The second point concerns the validity of two-dimensional solutions to predict significantly the main features of the flow patterns. Two-dimensional finite difference and spectral solutions are analysed with respect to the various results in a horizontal cylinder.

- The last point is related to the analysis of the fluid dynamics of physical vapour transport in a two-dimensional model with various axial wall temperature models. The analysis is made for a binary mixture in some horizontal and vertical situations.
<figure><img src="image_5.png" /><figcaption>Fig.4. Velocity fields in vertical midplane ($\phi=0°,180°$) for Ra$_T$=18,720, a=5 (L=10cm, R=1cm) and $\Omega=90°$: (a) experiments, (b) computations. From Smutek et al., 1985.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig.5. Three-dimensional velocity field (magnitude larger than 20% of the maximum) at various azimuths: (a) $\phi=45°,135°$ and (b) $\phi=90°,180°$; Ra$_T$=18,720, a=5. From Bontoux et al., 1986.</figcaption></figure>

### Three-Dimensional Convection in Tilted Cylinders

**Convection regimes in horizontal cylinders.** The core driven regime (CDR) is characterized by a parallel flow in the core and the w(r)-profile is S-shaped as predicted by relation (15a) (see Smutek et al., 1985, also in following Fig.18 for Ra$_T$=660). At Ra$_T$ ≈ 18,700 which is after the transition to the BLDR for a=5 (L=10cm, R=1cm, see Bontoux et al., 1986), the flow is inclined with respect to the axis (see Fig.4). The maximum w-velocity in the core has increased and shifted towards the walls (see following Fig.18 for Ra$_T$=18,720). As shown in Fig.5 secondary vortices develop near the end walls and boundary layers expand along the lateral walls (see Schiroky and Rosenberger, 1984, Smutek et al., 1985). In Fig.5 the fully three-dimensional velocity patterns are given at four regularly spaced azimuths : $\phi = \pi/4, 3\pi/4, \pi/2$ and $\pi$ and the velocity vectors are plotted only when they exceed some 10% of their maximum. In the regions hidden by the lateral (shaded) surfaces, they are only plotted at the last row in the r-direction. The circular section at the left of the cylinder corresponds to the cross-section at the plane nearest the cold end wall. The cross-flow structure superimposed on the main counterflows in the entire cavity is demonstrated in Fig.6 at various
<figure><img src="image_7.png" /></figure>

Fig.6. Flow structures in the (r,φ)-plane at various vertical cross-sections between the hot wall (z=2a) and the centre (z=a) for Ra$_T$=18,720 and a=5. From Smutek et al., 1983.

<figure><img src="image_8.png" /></figure>

Fig.7. Track of a particle released near the axis in the hot end region for Ra$_T$=18,720 and a=5.

<figure><img src="image_9.png" /></figure>

Fig.8. Schematic flow pattern in a horizontal cylinder in the BLDR. Primary, secondary flows. From Bontoux et al., 1986

cross-sections (r,φ) between the hot wall and the centre. The track of a fluid particule released at a point near the axis of the cylinder in the hot region is displayed in Fig.7. It reveals the occurence of both secondary vortices which set near the hot and cold end walls, and azimuthal transport in the axial direction (refer to the helicoidal flow structures emphasized in Fig.6). A sketch of the complete flow pattern is presented in Fig.8 at Ra$_T$≈18,700 (see Bontoux et al., 1986).

**Effect of the inclination Ω.** Experimental results by Rosenberger and co-workers are available in an a=5 cylinder (L=10cm, R=1cm) when the axis is tilted at an angle of 20°≤Ω≤150° with the gravity (see Schiroky and Rosenberger, 1984, also in Bontoux et al., 1986) at Ra$_T$=3580 and 18,720 which correspond to supercritical values in vertical cylinders heated from below (Ra$_{Tc}$≈431a when Ω=0° and a>>1). The computations were carried out for severe flow patterns at Ra$_T$=18,720 and 20°≤Ω≤90° (heating from below) with a 9x32x65-mesh.

The computed and experimental w-results in the core (z=5) are plotted vs. radius in Fig.9a. The comparison shows fairly good agreement (less than a 10% difference on the maximum of w when 30°≤Ω≤60°). Between Ω=90° and 30°, the w-profile changes from a Z-shape in BLDR at Ω=90° to a S-shape (again, as in CDR at Ω=90°) at Ω=30°. The
<figure><img src="image_10.png" /><figcaption>Fig.9. Effect of inclination Ω: (a) axial velocity (w) profiles vs radius at various Ω, computed (○, ▽, ◆) and experimental (●, ▼, ■) results; (b) maximal core w-velocity vs Ω, computed (△) and experimental (●, ■) results. From Hurford et al., 1982/84, Bontoux et al., 1986.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig.10. Velocity fields in the vertical midplane (φ=0°, 180°) at an inclination angle Ω=30° for Ra<sub>T</sub>=18,720 and a=5. Experiments (a,b) and computations (c) over 0≤z≤a (a) and 0≤z≤a/5 (b,c) from the cold end. (Refer to Fig.4.) From Hurford et al., 1982/84, Bontoux et al., 1986.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig.11. Three-dimensional velocity field (magnitude larger than 20% of the maximum) at various azimuths: (a) $\phi=90°, 180°$ and (b) $\phi=45°, 135°$; $Ra_T=18,720$, $a=5$, $\Omega=30°$. From Bontoux et al., 1986.</figcaption></figure>

variation vs. $\Omega$ of the maxima of $w$ is displayed in Fig.9a. Experimental results exhibit a dependence on both $\Omega$ and $Ra_T$. At $Ra_T=3580$ the variation of $w$ is well represented with $w_{max}(\Omega)/w_{max}(90°) \approx 1 - \sin(2\Omega - \pi)$ up to $\Omega=20°$. The experiments reveal a slight asymmetry at $\Omega \le 45°$ when $Ra_T=3580$ and $\Omega \le 60°$ when $Ra_T=18,720$. Also, the numerical results for $Ra_T=18,720$ are shown in Fig.9b to slightly underpredict the experimental variation.

The velocity fields in the vertical symmetry plane and in the vicinity of the cold end wall are shown in Fig.10a-c. The basic flow coming from the hot end wall separates into two parts in plane $\phi=0-\pi$ at the junction with the cold wall: - the main part is turned downward by the cold wall and at a certain distance from the endwall the flow becomes parallel again to the axis; - the second part is driven into a secondary counter-rotating vortex which develops in azimuth. The agreement is obvious between experiments (Fig.10b) and computations (Fig.10c) in both magnitude and direction; the discrepancies are small and mainly concern to the size and location of the vortices in the vertical plan.

The fully three-dimensional velocity field is given in Fig.11 at $\Omega=30°$ (same details as for $\Omega=90°$ in Fig.5). The velocity patterns show the occurrence of a weak flow everywhere in the plane ($\phi = \pi/2, 3\pi/2$) except near the end walls where boundary layers and vortices develop. The secondary co-rotating vortex which exists at the cold end for $\Omega=90°$, has disappeared at $\Omega=60°$. The counter-rotating vortex (see Fig.10b-c) is much smaller and expands in azimuth in a way which is similar to the
<figure><img src="image_13.png" /><figcaption>Fig.12. Isotherm patterns at Ω=30° and 60° for Ra<sub>T</sub>=18,720 and a=5 (isotherms spaced by ΔT/10). From Bontoux et al., 1986.</figcaption></figure>

vertical case (Smutek, 1984). Except in the end regions the flow is mainly parallel to the axis as shown in Fig.11 for $\phi = \pi / 4$, $3\pi / 4$ and $\pi$.

The isotherm patterns displayed in Fig.12 exhibit distortions close to the end walls and in the core which result from the increase of transport at Ω<90° (see following Fig. 19 for Ω=90°). The isotherms also reveal the onset of the secondary counter-rotating vortices shown in Fig.10 at Ω<60°.

**Axisymmetric and asymmetric regimes in vertical cylinders (Ω=0°).** The first transition from rest to steady convection in a vertical cylinder heated from below is predicted by linear stability analyses (see Charlson and Sani, 1971, Gershuni and Zhukhovitski, 1976, Buell and Catton, 1983) to be dominated by the axisymmetric mode in flat cylinder (a<a<sub>c</sub>) and by the antisymmetric mode in long cylinder (a>a<sub>c</sub>). The study of these regimes is detailed elsewhere (see Smutek, 1984, Bontoux et al., 1986), however, some points are reported here.

Axisymmetric flows were computed at small a=0.5 and at a=0.625≈a<sub>c</sub>. For this last value of aspect ratio a, both axisymmetric and asymmetric patterns were obtained depending on the initial disturbances imposed to the fluid at rest. Also, as shown in Fig.13, the axisymmetric solutions can involve either an upstream or a downstream at the centre.

Attention was focused on the asymmetric solution in an a=1 cylinder above critical Ra<sub>Tc</sub>. At Ra<sub>T</sub>/Ra<sub>Tc</sub>≈5.3 secondary counter-rotating vortices are superimposed on the basic flow in the symmetry plane ($\phi = 0, \pi$) and develop in azimuth into a four-eddy pattern in the plane ($\phi = \pi / 2, 3\pi / 2$) as shown in the two-dimensional velocity fields displayed in Fig.14 at Ra<sub>T</sub>=6250. The three-dimensional velocity field is displayed in Fig.15 which emphasizes the azimuthal expansion of the plane rolls shown in Fig.14. A schematic pattern is proposed in Fig.16. The occurrence of such patterns is confirmed by solutions obtained with insulated side walls (see Müller et al., 1984).
<figure><img src="image_14.png" /><figcaption>Fig.13. Axisymmetric velocity field at a=0.625, Ω=0° and Ra<sub>T</sub>=2000: (a) ascending and (b) descending flow in the centre. From Smutek, 1984, and Bontoux et al., 1986.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig.14. Asymmetric solution at a=1, Ω=0° and Ra<sub>T</sub>=6250. Velocity fields in various planes (r,z). From Smutek, 1984.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig.15. Three-dimensional velocity field (magnitude larger than 20% of the maximum) at various azimuths: (a) φ=45°,135° and (b) φ=90°,180°; Ra<sub>T</sub>=6 a=1, Ω=0 From Bontoux et al., 1986.</figcaption></figure>

<caption>Fig.16. Schematic flow pattern for the asymmetric regime in a vertical cylinder heated from below. From Bontoux et al., 1986.</caption>
Analytical solutions for the maximum values of the horizontal w-
velocity component can be derived from (14) and (15a) respectively as

$$
w^{2D}_{\max} = k_1 Ra_T / 9\sqrt{3} \quad (26)
$$

and

$$
w^{3D}_{\max} = k_1 Ra_T / 12\sqrt{3} \quad (27)
$$

These solutions are plotted for a=5 with the conduction condition
k₁=(2a)⁻¹ in Fig.17 together with various experimental and numerical
results. The numerical results are obtained with the three-dimensional
(3-D) finite difference solutions (FD3D), the two-dimensional (2-D)
hermitian finite difference (HFD2D) and Tau-Chebyshev spectral (TC2D)
solutions. In addition, Fig.17 contains plots of (27) with evaluations of k₁
obtained for the BLDR and also with the third-order expansion term in
(Ra_Tk₁) proposed by Schiroky and Rosenberger (1984) (see Bontoux et
al., 1986). The analysis of these curves results in the following
observations:

(i) The experimental and 3-D numerical results agree well over the
range covered by the numerical data (74<Ra<30,000). The 2-D results
show significant deviations from these 3-D results not only, as expected,
in the CDR but also in the BLDR.

(ii) At low Ra<sub>T</sub> (in the CDR) the 3-D analytical solution (27) yields good predictions. The 2-D solution (26), though correctly reflecting the linear dependence of w<sub>max</sub> on Ra<sub>T</sub>, overestimates by 1/3 with respect to the cylinder.

(iii) At high Ra<sub>T</sub> (i.e. in the BLDR) the experimental and 3-D numerical results for the core velocity vary as Ra<sub>T</sub><sup>1/2</sup>, parallel to Gill (1966)'s relation for boundary-layer. The velocity in the cylinder is about 1.5 times larger than the 2-D numerical solution, which however shows a realistic Ra<sub>T</sub>-dependence. The velocities obtained from the analytical solutions derived by Schiroky (1982) (see Bontoux et al., 1986) are too low and their Ra<sub>T</sub>-dependence is nearly Ra<sub>T</sub><sup>0.4</sup>.

(iv) If we define a critical Ra<sub>Tc</sub> for the transition between the CDR and BLDR, then we find a factor of about 3 for a=5, between the 2-D (Ra<sub>Tc</sub><sup>2D</sup>~2000) and 3-D (Ra<sub>Tc</sub><sup>3D</sup>~6000) behaviours.

The above points are further illustrated by w(r)-profiles for two
Ra values in the CDR and BLDR, and by streamline and isotherm patterns
compared to 3-D velocity fields and isotherm patterns. In Fig.18a one
sees that in the CDR (Ra<sub>T</sub>=660) the 3-D analytical results agree well with
the experimental and numerical values. The 2-D numerical values lie too
<figure><img src="image_17.png" /><figcaption>Fig.17. Maximum value of dimensionless axial velocity (core velocity) vs Ra$_T$ for $\Omega=90^\circ$: comparison of experimental, analytical and numerical results.</figcaption></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>Fig.18. Core velocity (w) profiles for: (a) Ra$_T$=660 and (b) Ra$_T$=18,720. $\Omega=90^\circ$, a=5.
From Bontoux et al., 1986.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig.19. Computed (FD3D) velocity fields and isotherms in the vertical midplane of a cylinder. Computed (HFD2D and TC2D) streamlines and isotherms in a cavity. a=5 and Ω=90°. Comparison for various Ra$_T$: (a) Ra$_T$=3580, (b) Ra$_T$=8860 and (c) Ra$_T$=18,720.</figcaption></figure>

high. At Ra$_T$=18,720 (Fig.18b), the 2-D numerical results and 3-D analytical velocity profiles differ strongly from each other and from experiments. There are also distinctions into "S-shapes" and "Z-shapes". The evolution of the Z-shape from the S-shape occurs slower for 2-D than for 3-D solutions. One should point out the good agreement between experiment and 3-D numerical solutions.

The streamline and isotherm patterns are presented with the velocity fields computed in the vertical plane of the cylinder in Fig.19 near and above the transition to the BLDR. The agreement is rather good
<figure><img src="image_21.png" /><figcaption>Fig.20. Three-dimensional velocity field (magnitude larger than 20% of the maximum) over a distance of a/5 from the hot wall. Ra$_T$=18,720, a=5. From Bontoux et al., 1986.</figcaption></figure>

between the (21x41) finite difference solutions and the (12x27) Tau-Chebyshev solutions. The streamlines in the cavity remain nearly parallel to the horizontal walls up to Ra$_T$=18,720 while the velocity fields are already inclined at Ra$_T$=8860 in the cylinder. Also, the onset of secondary vortices near the vertical walls is revealed in both 2-D and 3-D computations but they arise differently near the cold wall : below the axis in the cavity, and above in the vertical plane of the cylinder. The distortion of the isotherms near the end wall is more important in the cavity than in the cylinder. This reflects the three-dimensional transport. As suggested by the experimental results of Schiroky and Rosenberger (1984) and observed in the 3-D numerical simulation by Smutek et al. (1985) and Bontoux et al. (1986), the U-turn of the flow occurs there over more than 180° in azimuth as displayed in Fig.20.

## Simulation of Physical Vapour Transport in Cavities

The numerical solution of the governing system (4-7) is considered for PVT with both spectral Tau-Chebyshev method (TC2D) and mixed $O(h^2, h^4)$ finite-difference method (HFD2D) detailed in previous sections. Chebyshev-polynomial expansions are well adapted when the gradient zones are located near the boundaries. They are used with the 1-G thermal model (Fig.2a) and the linear Fick's law (8) at the interfaces. For a=4 the computations are carried out with no more than 17x17 Chebyshev components when $Ra_T \le 6.10^3$ and $|Ra_S| \le 2.10^3$ (Elie et al., 1986). For the multizone thermal models, 2-Z and 3-Z, the mixed finite differences are used with a 21x81 uniform mesh (Extrémet et al., 1986). Also, the simple analytical model (10) is taken at the interfaces.

**Fick's Law and 1-G Thermal Model.** Negative solutal Rayleigh numbers correspond to binary mixtures where the heavier component ($M_A \ge M_B$) is the active one (A) in PVT process ($S_{A2} \ge S_{A1}$). Then, solutal and thermal buoyancy terms act in opposite directions giving rise to more complex flow patterns than with $Ra_S > 0$. In Fick's model (8) the values for E (and F) are derived from relation (9) with $1<S_{A0}<0.5$ and $\Delta S_A < 0.5$. The computations are carried out for steady solutions in cavities with $1 \le a \le 4$ at horizontal ($\Omega=90°$) and vertical inclinations (with heating from below, $\Omega=0°$, and from above, $\Omega=180°$).
<figure><img src="image_22.png" /><figcaption>Fig.21. Streamlines for various E at Ra$_{S}$=-12.5, Ra$_{T}$≈0, a=4. (Elie, 1984)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig.22. Streamline, isotherm, iso-concentration patterns and velocity fields (A-B components) for E=10 at Ra$_{S}$=-12.5, Ra$_{T}$≈0, a=4.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig.23. Streamline, isotherm, iso-concentration patterns and velocity fields (A-B components) for E=1.05 at Ra$_{S}$=-12.5, Ra$_{T}$≈0, a=4. From Elie, 1984, also Elie et al., 1986.</figcaption></figure>

(i)Horizontal cavities ($\Omega$ =90°). In horizontal cavities the temperature and concentration gradients between the source and the sink are perpendicular to the gravity vector. As above mentioned a basic flow, then, exists even at very low Ra and can be descending at the (hot) source depending on the density gradient due to dominant temperature or concentration change. Also the mass flux of component A superimposes to thermal and solutal convection an advective diffusive motion through the cavity. A variety of such complex flows is studied with Pr=Le=1 at a=4.
<figure><img src="image_25.png" /><figcaption>Fig.24. Computed concentration profiles along the axis for various E (related to a Peclet number) at Ra$_S$=-12.5, Ra$_T$≈0, a=4. From Elie, 1984. (Refer also to Greenwell et al., 1981.)</figcaption></figure>

When $M_A > M_B$ the streamline pattern (at Ra$_S$=-12.5 and Ra$_T$≈0 in Fig.21) is strongly affected by an increase of the interfacial mass flux at the source (as E varies from 10 to 1.05). The basic cell is descending at the source and ascending at the sink. The mass flux of A organizes into a forced flow which develops from the source to the sink below the solutal cell. When the flux rate increases, the cell is progressively damped towards the sink. Typical velocity fields are given in Fig.22 and 23 with streamline, isotherm and iso-concentration patterns, at small (E=10) and large (E=1.05) flux rates and for F=10. When E=10 (Fig.22) the motion of B is mainly driven by the solutal cell, while component A is dominantly advected in the lower part of the cavity by the net flow between the interfaces. When E=1.05 (Fig.23) the advection from the source dominates. The concentration varies then strongly near the sink (see the concentration profiles at various E in Fig.24). The mass average and A-component velocity fields are very similar and the motion of B occurs in the regions where A-velocities are the lowest, i.e. in the solutal cell and near the horizontal walls.

When $M_A \approx M_B$ (Ra$_S$≈0) the basic cell is driven by the temperature changes. The basic cell is ascending at the hot source and descending at the cold sink (see Fig.25). Then, the net flow develops in the upper part of the cavity and gets more and more confined as Ra$_T$ is increased. At Ra$_T$=5860 the streamlines suggest the onset of a secondary cell near the sink (crystal). (See in Fig.19 for a mono-component gas). Informative patterns and velocity fields are presented in Fig.26 to 28 for Ra$_T$=19.5, 586 and 1950. The flow of A is mainly driven by the net flow. The transport of B is associated with the thermal cell. Strong gradients develop at Ra$_T$=1950 (Fig.28) near the vertical walls.
<figure><img src="image_26.png" /><figcaption>Fig.25. Streamline patterns when Ra$_T$ varies at Ra$_S$≈0, E=5.27, a=4.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig.26. Streamline, isotherm, iso-concentration patterns and velocity fields (A-B components) for E=5.27 at Ra$_T$=19.5, Ra$_S$≈0, a=4.</figcaption></figure>

When Ra$_T$ ≈ |Ra$_S$|/Le there is a direct competition between the thermal and solutal buoyancy forces. At a large flux rate (E=1.20) the results are displayed in Fig.29 for Ra$_T$=11.7 and Ra$_S$=-12.5. The net flow dominates but the opposite effects of concentration and temperature gradients are still revealed by slight distortions in the streamline pattern. In the following section the competition will be emphasized for a multizone thermal model.

(ii) Vertical cavities heated from below (Ω=0°). In a cylindrical enclosure filled with monocomponent gases, a steady antisymmetric instability is predicted to onset at a=1 by both stability analyses and computations (see above sections). Computed 2-D solutions involving a net flow (E=2.3) are found, however, to be axisymmetric at Ra$_T$=3800 and Ra$_S$≈0 with Pr=1 and Le=0.114 (Fig.30) in agreement with previous results by Markham et al. (1980, 1981). The flow is driven upward along the axis and two cells develop on both sides due to temperature gradients.
<figure><img src="image_28.png" /><figcaption>Fig.27. Streamline, isotherm, iso-concentration patterns for E=5.27 at Ra$_T$=586, Ra$_S$≈0, a=4.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig.28. Streamline, isotherm, iso-concentration patterns for E=5.27 at Ra$_T$=1950, Ra$_S$≈0, a=4.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig.29. Streamline, isotherm, iso-concentration patterns for E=1.20 at Ra$_T$=11.7, Ra$_S$=-12.5, a=4.</figcaption></figure>

(iii) Vertical cavities heated from above (Ω = 180°). The computations are made again for a=1 with Ra$_T$=500, Ra$_S$=-1660, Pr=1, Le=0.18, E=2.3 (S$_A$0=0.99 and Δ S$_A$=0.009). Due to dominant advection diffusion an axisymmetric downward motion develops along the axis (Fig.31) with two small side cells near the vertical walls and driven by the density gradients. When a=4 with Ra$_T$=541 and Ra$_S$=-1475 with Pr=1, Le=0.18 and same E=2.3 (Fig.32) the flow pattern is noticeably different. The fluid does not flow longer downward along the axis but now along the side walls with two recirculating cells in between symmetrically to the axis.
<figure><img src="image_31.png" /><figcaption>Fig.30. Streamline, isotherm, iso-concentration patterns for E=2.3 at Ra$_T$=3800, Ra$_S$≈0, Le=0.114, Pr=1, a=1 and Ω=0°.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Fig.31. Streamline, isotherm, iso-concentration patterns for E=2.3 at Ra$_T$=500, Ra$_S$=-1660, Le=0.18, Pr=1, a=1 and Ω=180°.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Fig.32. Streamline, isotherm, iso-concentration patterns for E=2.3 at Ra$_T$=541, Ra$_S$=-1475, Le=0.18, Pr=1, a=4 and Ω=180°.</figcaption></figure>
<figure><img src="image_34.png" /></figure>

Fig.33. Variation of the maximum axial velocity with $a\gamma$ in the conduction regime for various computed values of a and the analytical laws (26) and (17). From Extrémet et al., 1986.

**Multizone Thermal Models and Horizontal Cavities (Ω=0°).** The analytical solution (26) with $k_1$ derived from (17a-b) have been compared in Fig.33 to the maximum velocities computed for 2-Z thermal models (Fig.2b) at various a and $a\gamma$. The agreement is fairly good as the computed solutions shows a ($a\gamma$)-dependence when $a\gamma>1$ and are independent on $a\gamma$ when $a\gamma<1$. When $a\gamma<1$ the computed values, however, stay substantially below the prediction by about 60% at most (refer to discrepancy factor $\beta$ in (18)). Also, the analytical solution (19) is shown elsewhere to predict correctly the order of magnitude (at least) of the velocities for 3-Z thermal models (Fig.2c) and binary mixtures (see Extrémet et al., 1986).

With a 3-Z thermal model more than with a 1-G model, completely different flow patterns can be expected if $Ra_T \gg Ra_S$/$Le$ (thermal convection dominated regime), $Ra_T \ll Ra_S$/$Le$ (solutal convection dominated regime) or $Ra_T \approx Ra_S$/$Le$ (mixed regime). For the mixed regime the buoyancy terms have coefficients of comparable magnitudes everywhere and the density gradients determines locally the resulting ascending or descending motion. The flow patterns also depend on the strength of interfacial fluxes modelized with relation (10) and characterized by $v_d$.

(i) Without interfacial flux ($v_d=0$) for $a=4$ and $\Theta_m = (T_m - T_1)/\Delta T=0.75$, the flow in the mixed regime ($Ra_S \approx 40$, $Ra_T \approx 19.5$, $Le \approx 1.4$) is composed of two counter-rotating cells. (See Fig.34 and compare with the solution for $Ra_S \approx 0$). The larger one, extending over 5/6 of the cavity, is driven by dominant solutal buoyancy. The motion is clockwise and the effect of the small temperature gradient (of width $\delta_2$) which drives opposite buoyancy forces, is revealed by the streamline pattern.
<figure><img src="image_35.png" /><figcaption>Fig.34.Streamlines, isotherms, iso-concentrations and axial temperature profiles for a=4, Ra<sub>T</sub>=19.5, Le=1.4, and: (a) Ra<sub>S</sub>=0, (b) Ra<sub>S</sub>=-40. From Extrémet et al., 1986.</figcaption></figure>

<caption>Fig.35.Streamlines, iso-concentrations for a=4, Ra<sub>T</sub>=19.5, Le≈1.4 and lv<sub>d</sub>=2 at various Ra<sub>S</sub>: (a) Ra<sub>S</sub>=-40, (b) Ra<sub>S</sub>=-21, (c) Ra<sub>S</sub>=-3. From Extrémet et al., 1986.</caption>

(ii)With an interfacial flux (ν<sub>d</sub> ≠ 0) the modifications of the streamline pattern are displayed in Fig.35. The forced flow entering the cavity from the source, develops under the "solutal cell" then above the "thermal cell" (Fig.35a). In Fig.35b-c we show the flow patterns when Ra<sub>S</sub> decreases. The forced flow mostly runs in the lower part of the domain when Ra<sub>T</sub> < |Ra<sub>S</sub>|/Le at Ra<sub>S</sub> ≈ -40 (except in the close vicinity of the sink, see Fig.35a). When Ra<sub>T</sub> > |Ra<sub>S</sub>|/Le (Ra<sub>S</sub> ≈ -3) it moves entirely in the upper part of the cavity (Fig.35c). In this case, the secondary cells in the bottom of the cavity are anti-clockwise and mostly driven by the thermal buoyancy force in the two gradient-zones. In Fig.35b, the thermal and solutal buoyancy balances as Ra<sub>T</sub> ≈ |Ra<sub>S</sub>|/Le (Ra<sub>S</sub> ≈ -21). Solutal buoyancy still dominates near the source but the basic "solutal" cell is, then, cut into two co-rotating cells by the second thermal gradient (ϑ<sub>2</sub>). Near the sink, the "thermal" cell corresponding to the first gradient (ϑ<sub>1</sub>) extends over the entire height. The forced flow visualized by one streamline, is shown to slide between these three cells.
CONCLUSIONS

The paper is an amalgam of various (published and un-published)
numerical works which common purpose is the understanding of the
fluid dynamics and the heat and mass transfer for vapour crystal
growth applications in enclosures. The experimental studies carried out
in the group of Prof. F. Rosenberger were used for the validation of the
three-dimensional numerical solutions in cylinders. The results concern
the following points:

(i) The complex flow patterns were analysed for various inclinations with respect to the gravity vector. The transition from conduction dominated regime (core driven regime, CDR) to boundary layer driven regime (BLDR) is characterized for a horizontal cylinder. A good agreement is obtained between computations and experiments in long tilted cylinders for the BLDR, and between computations and stability analyses at the onset of convection in vertical cylinders (at moderate aspect ratios).

(ii) The limitations of the 2-D model to predict the main features of the flow were studied for a mono-component gas and aspect ratio $a=5$. The velocity scales can be under- or over-estimated by 30 to 50% depending on $Ra_T$. The $Ra_T$ for the transition between CDR and BLDR for the core velocity is underestimated with the 2-D solutions by a factor of 3 ($Ra_{Tc}^{2D} \approx 2000$, $Ra_{Tc}^{3D} \approx 6000$).

(iii) The analysis of a vapour growth process (PVT) was made with 2-D solutions and using mass-flux models at the interfaces. Both analytical approximations and direct simulations were considered for one-gradient and multizone (wall) temperature profiles. Very complex flow patterns were emphasized for various $Ra_T$ and $Ra_S$. Agreement was obtained also with other numerical results obtained with different methods.

Experiments are in process in the laboratory of F. Rosenberger and data should be soon available for comparison with numerical solutions. Some effects connected with the presence of a seed located (for cristallization) at the side walls of a horizontal ampoule were discussed by Extrémet et al. (1986). Also, 3-D solutions were computed for multizone devices (Extrémet, 1986). A recent modelling has been used in the interface conditions which introduce a coupling between convection in the fluid and surface kinetics (see Zappoli and Elie, 1984, and Zappoli, 1986).

**Acknowledgments:** The authors want to gratefully acknowledge Dr. S.S. Leong for use of code CYL3D and Dr. J.M. Lacroix for participation in the elaboration of graphics softwares. They wish to thank Prof. R. Cadoret and Dr. J.C. Launay for focusing their attention on the fluid dynamics problems involved in their growth experiments. Collaboration with Dr. B. Zappoli is also acknowledged. Financial supports by the D.R.E.T. (G.6), the C.N.E.S. (Division "Matériau") and the C.N.R.S. (GRECO-64). Computing centers: C.C.V.R. (CRAY1/S), C.N.U.S.C. and P.A.-Saint Charles (K. Dang Quoc, G. Barisone, D. Pardini, J.L. Zubillaga).
REFERENCES

Adam Y., 1977, J. Comput. Phys., 24, 10-22.

Bejan A., Tien C.L., 1978, Int. J. Heat Mass Transfer, 21, 701.

Bejan A., Tien C.L., 1978, J. Heat Transfer, 100, 641-647.

Birikh R.V., 1966, J. Appl. Math. & Tech. Phys., 30, 432.

Bondet de la Bernardie B., 1980, Doct. Thesis, Univ. Aix-Marseille III.

Bontoux P., Forestier B., Roux B., 1978, J. Méc. Appl., 2, 3, 291-316.

Bontoux P., Bondet de la Bernardie B., Roux B., 1981, in Numer. Meth. Coupled Prob., Ed. Hinton, Bettess, Lewis, Pineridge, 1018-1030.

Bontoux P., Smutek C., Roux B., Lacroix J.M., 1986, J. Fluid Mech., 169, 211-227.

Bontoux P., Roux B., Schiroky G.H., Markham B.L., Rosenberger F., 1986, Int. J. Heat Mass Transfer, 29, 227-240.

Bontoux P., Smutek C., Roux B., Extrémet G.P., Schiroky G.H., Hurford A.C., Rosenberger F., 1986, in Numer. Meth. Non-Linear Prob., 3, Ed. Taylor, Owen, Hinton, Damjanic, Pineridge, 1102-1115.

Bontoux P., Smutek C., Randriamampianina A., Roux B., Extrémet G.P., Hurford A.C., Rosenberger F., 1986, in Advances in Space Research, Ed. Legros, Pergamon.

Bontoux P., Smutek C., Roux B., Hardin G., Sani R., 1986, Three-Dimensional Buoyancy Driven Flows in Cylindrical Cavities with Differentially Heated End Walls. Part II : Vertical Cylinders (submitted to J. Fluid Mech.).

Brisson P., 1981, Ing.-Doct. Thesis, Univ. Clermont 2, France.

Buell J.C., Catton I., 1983, J. Heat Transfer, 105, 2, 255-260.

Carruthers J.R., 1977, in Prep. and Prop. of Solid State Materials, 3, Ed. Wilcox and Lefever, Dekker.

Charlson G.S., Sani R., 1971, Int. J. Heat Mass Transfer, 14, 2157-2160.

Cormack D.E., Leal L.G., Imberger J., 1974, J. Fluid Mech., 65, 209.

Cormack D.E., Leal L.G., Seinfeld J.H., 1974, J. Fluid Mech., 65, 231.

Cooley J., Tukey J.W., 1965, Math. Comp., 19, 297.

De Vahl Davis G., 1979, Numer. Heat Transfer, 2, 261.

De Vahl Davis G., 1986, Proc. 8th Int. Heat Transfer Conf., ASME, 101-109.

Elie F., Chikhaoui A., Randriamampianina A., Bontoux P., Roux B., 1983, Notes Numer. Fluid Mech., 7, Vieweg, 57.

Elie F., 1984, Doct. Thesis, Univ. Aix-Marseille II, France.

Elie F., Randriamampianina A., Bontoux P., Extrémet G.P., Roux B., 1986, in Numer. Meth. Non-Linear Prob., 3, Ed. Taylor, Owen, Hinton, Damjanic, Pineridge, 738-753.

Extrémet G.P., Bontoux P., Roux B., 1986, Int. J. Heat Fluid Flow (in print)

Extrémet G.P., Roux B., Bontoux P., Elie F., 1986, J. Crystal Growth (in print)

Extrémet G.P., Bontoux P., Roux B., 1986, in Advances in Space Research, Ed. Legros, Pergamon.

Extrémet G.P., 1986, "Numerical Modellings for Multizone Physical Vapour Transport. Two- and Three-Dimensional Solutions", to be published in Materials Sciences under Micro-Gravity, ESA Scien. & Tech. Publications.
Gershuni G.Z., Zhukhovitskii E.M., 1976, "Convective Stability of Incompressible Fluids", Keters/Wiley.

Gill A.E., 1966, J. Fluid Mech., 26, 3, 515-536.

Gottlieb D., Orszag S.A., 1977, CBMS-NSF Conf. Series Applied Math.

Greenwell D.W., Markham B.L., Rosenberger F., 1981, J. Crystal Growth, 51, 413.

Haidvogel D.B., Zang T., 1979, J. Comp. Phys., 30, 167.

Hart J. E., 1972, J. Atmos. Sci., 29, 687-697.

Hart J. E., 1983, Int. J. Heat and Mass Transfer, 26, 1069-1074.

Hindmarsh A.C., LLNL Rep., UCID-30130, 1976 & UCRL-88007, 1982

Hirsh R.S., 1975, J. Comp. Phys., 10, 90-109.

Hurford A.C., Schiroky G.H., Rosenberger F., 1982/84-Reports, Univ. Utah, Salt Lake City, USA.

Imberger J., 1974, J. Fluid Mech., 65, 247.

Jhaveri B.S., Markham B.L., Rosenberger F., 1981, Chem. Eng. Comm., 13, 65.

Jhaveri B.S., Rosenberger F., 1982, J. Crystal Growth, 57, 57.

Joseph D.D., 1976, "Stability of fluid motions II", Springer Tracts in Natural Philosophy, 28.

Klosse K., Ullersma P., 1973, J. Crystal Growth, 18, 167-174.

Krause E., 1972, V.K.I. Lecture Series on Numer. Meth. Fluid Dynamics, Rhode-St-Genèse, Belgium.

Krause E., 1974, Numerical solutions of the Navier-Stokes equations", L.N. I.C. for Mech. Sc., Udine, Italy.

Krause E., Hirschel E.H., Kordulla W., 1976, Computers and Fluids, 4, 2, 77-92.

Launay J.C., 1982, J. Crystal Growth, 60, 185-190

Le Bail R.C., 1972, J. Comput. Phys., 9, 440,

Leong S.S., De Vahl Davis G., 1979, Numer. Meth. Thermal Prob., Ed. Lewis & Morgan, Pineridge, 287

Leong S.S., 1983, PhD thesis, Univ. of NSW, Kensington, Australia.

Lhomme B., Morgenstern J., Quandalle P., 1982, Euromech 159, Nice.

Loc T.P., Daube O., 1978, J. Mécanique, 17, 5, 651-678.

Mallinson G.D., De Vahl Davis G., 1973, J. Comput. Phys., 12, 436.

Markham B.L., Rosenberger F., 1980, Chem. Eng. Comm., 5, 287.

Markham B.L., Greenwell D.W., Rosenberger F., 1981, J. Crystal Growth, 51, 426.

Markham B.L., Rosenberger F., 1984, J. Crystal Growth, 67, 241.

Müller G., Neumann G., Weber W., 1984, J. Crystal Growth, 70, 78.

Ouazzani J., Peyret R., 1983, Notes Numer. Fluid Mech., 7, Vieweg, 57.

Ouazzani J., 1984, Doct. Thesis, Univ. Nice, France.

Olson J.M., Rosenberger F., 1979, J. Fluid Mech., 92, 4, 631-642.

Omaly J., Robert M., Cadoret R., 1981, Mat. Res. Bull., 16, 1261.

Omaly J., Robert M., Brisson P., Cadoret R., 1983, Nuclear Inst. and Meth., North-Holland, 213, 19-26.

Orszag S.A., 1971, J. Fluid Mech., 149, 75.

Ostrach S., Loca R.R., Kumar A., 1980, Natural Convection in Enclosures, Ed. Catton & Torrance, ASME, HTD-8, 1.
Ostrach S.,1982, Ann. Rev. Fluid Mech., 14, 313-345.

Pamplin B., 1980, Int. Series Sciences Solid State, 16, Ed. Pamplin, Pergamon.

Peaceman D.W., Rachford H.H., 1955, J. Soc. Ind. Appl. Math., 3, 28.

Petzold L.R.,1980, Sandia Nat. Lab. Rep. SAND80-8230.

Peyret R., Taylor T.D., 1982, "Computational Methods for Fluid Flow", Springer-Verlag.

Pimputkar M., Ostrach S., 1981, J. Crystal Growth, 55, 614-646.

Randriamampianina A.,1984, Doct. Thesis, Univ. Aix-Marseille II, France.

Randriamampianina A., Bontoux P., Roux B., Argoul P., 1985, Notes Numer. Fluid Mech., Vieweg, 13, 302.

Rosenberger F.,1979, Phys. Chem. Hydrodyn., 1, 1.

Rosenberger F.,1979, Fundamentals of Crystal Growth I, Solid State Sci., 5, Springer.

Rosenberger F., 1982, Convective Transport and Instability Phenomena, Ed. Zierep & Oertel, Verlag-Braun, 469.

Roux B., Grondin J.C., Bontoux P., Gilly B., 1978, Numer. Heat Transfer, 1, 331-349.

Roux B., Bontoux P., Ta Phuoc Loc, Daube O., 1980, in L.N. in Math. 771, Springer, 450-468.

Samuels M.R., Churchill S.W.,1967, A.I.Ch.E. J., 13, 77.

Shih T.S.,1981, Int. J. Heat Mass Transf., 24, 1295-1303.

Schiroky G.H., 1982, Ph.D. Thesis, Univ. Utah, Salt Lake City, USA.

Schiroky G.H., Rosenberger, F.,1984., Int. J. Heat Mass Transf., 27, 587

Smutek C., Bontoux P., Roux B., Schiroky G.H., Hurford A., Rosenberger F.,1985, Numer. Heat Transfer, 8, 613,

Smutek C.,1984, Doct. thesis, Univ. Nice, France.

Smutek C., Roux B., Bontoux P., De Vahl Davis G.,1983, Notes Numer. Fluid Mech., Vieweg, 7, 338,

Solan A., Ostrach S.,1979, Preparation and Properties of Solid State Materials, 4, 63-110.

Temperton C.,1983, J. Comp. Phys., 52, 81.

Tichy J., Gadgil A., 1982, J. Heat Transfer, 104, 103.

Zappoli B., Elie F.,1984, "On Boundary Conditions for Hydrodynamic Equations in Reduced Gravity", 35th Int. Astronautical Federation Conf., Lausanne, Switzerland.

Zappoli B., 1986, "On the Interaction between Convection and Surface Reactions in Rectangular Horizontal Enclosures", (submitted to J. Crystal Growth).