- Email: [email protected]
Modelling of transport phenomena in a low-pressure CVD reactor
ARTICLE IN PRESS
Journal of Crystal Growth 267 (2004) 598–612
Modelling of transport phenomena in a low-pressure CVD reactor Arnab K. Dea, K. Muralidhara,*, V. Eswarana, V.K Wadhawanb a
Department of Mechanical Engineering, Indian Institute of Technology, IIT Kanpur, Kanpur (UP) 208 016, India b Laser Materials Division, Center for Advanced Technology, Indore 452 013, India Received 13 April 2004; accepted 16 April 2004
Communicated by J.J. Derby
Abstract The present work reports the analysis of transport phenomena in a low-pressure chemical vapour deposition (CVD) reactor. A two-dimensional unsteady model has been used to compute flow, chemical reactions and mass transfer rates within the reactor. The formation of zinc sulphide from the reaction of zinc vapour and gaseous hydrogen sulphide has been analysed. Two geometric designs of the reactor have been studied for assessing the overall deposition characteristics. The mathematical model comprises of the conservation equations of mass, momentum and the mass fractions of the reactant and product species. The governing partial differential equations have been solved by the finite volume method. The deposition surface has been taken to be chemically passive. There is an uncertainty in the boundary condition for the depositing species on the passive surface. The difficulty has been resolved by utilizing a homogenous Dirichlet condition at the surface. The present approach yields an upper bound on the deposition rate. The numerical results obtained show that the species concentration field, correlates closely to the flow field in the reactor. Over the range of Reynolds numbers studied namely 10–200, a recirculation pattern adjacent to the side walls is revealed and a stagnation zone forms ahead of the block. Reactant species mix predominantly in the zone close to the inflow plane and the product species are transported by advection. The reactor with a convex substrate streamlines the flow better as compared to the concave. As a result, higher concentrations of the product species reach the substrate. At a Reynolds number of 100, the deposition rate for the convex substrate was seen to increase by a factor of two with respect to the concave. r 2004 Elsevier B.V. All rights reserved. PACS: 81.10.Bk; 81.15.Gh; 02.60.Cb; 02.70.Fj; 51.10.ty; 51.20.td Keywords: A1. Computer simulation; A1. Convection; A1. Diffusion; A1. Fluid flows; A1. Mass transfer; A2. Growth from vapor
*Corresponding author. Tel.: +91-5122597182; fax: +91-5122597408. E-mail addresses: [email protected] (A. K. De), [email protected] (K. Muralidhar), [email protected] (V. Eswaran), [email protected] (V.K. Wadhawan). 0022-0248/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2004.04.036
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
599
Nomenclature Cf Di D Mi Mc Pem Re rb kf kr s Uc Vr x; y; z
Zb
skin friction coefficient, tw =12rc Uc2 binary diffusion coefficient of species i in argon, m2 =s diameter of the reactor molecular weight of the ith species molecular weight of the carrier gas, namely argon Peclet number associated with mass transfer, Uc D=Di reynolds number, rUc D=m dimensionless radius of the substrate forward reaction rate, m3 /mol s backward reaction rate, m3 /mol s thickness of the deposited zinc sulphide, mm velocity of the central jet, m/s central-to-peripheral jet velocity ratio Cartesian coordinates in the physical domain with x pointing in the flow direction dimensionless distance of the block from the inflow plane
1. Introduction Chemical vapour deposition (CVD) is an industrial process for growing thin layers (0.01– 10 mm) of semiconductor materials, such as epitaxial silicon and gallium arsenide. Thin films are technologically important, especially in the field of microelectronics and as protective coatings for optical elements. As against thin films, the application of CVD to the bulk growth of monolithic compounds such as zinc sulphide and cadmium telluride has been considered in the present work. Infrared sensors are an important application where such materials are required. The deposited material has thicknesses in the range 1–10 mm. The reactor configuration specifically analysed operates in the batch-mode, and employs low-gas pressures and high temperatures. The CVD reactor is a large tube with inflow on one side and a deposition surface on the other. It contains nozzles at the inlet plane for separate
Greek letters a injection angle of the peripheral jet m dynamic viscosity of argon, Pa s stoichiometric coefficient for ith species ni in the chemical reaction oi mass fraction of ith species oA;o mass fraction of hydrogen sulphide at the inlet Rfj forward rate of reaction of jth species, mol m3 s r Rj backward rate of reaction of jth species, mol m3 s rc density of the carrier gas, namely argon, kg/m3 rZnS density of zinc sulphide, kg/m3 tw wall shear stress, N/m2 Subscript o inflow plane Superscripts g gas phase s surface properties
injection of the two reactants, each mixed with an inert carrier gas. The substrate is a flat or a curved surface that is placed at normal incidence to the incoming flow. The flow bypass occurs at the periphery of the substrate. The incoming gases chemically react, resulting in fine particles of the product. These are carried by the main flow and deposited on the substrate surface. The particles, in addition may be adsorbed and later mobilized by diffusion within the substrate. The other gaseous products are exhausted along with the main flow from the reactor. In the present study, the reactants are zinc vapour and hydrogen sulphide, while the carrier gas is argon. The product species of importance to the present work is zinc sulphide. Deposition of zinc sulphide can occur at all the solid surfaces of the reactor. The deposition on the substrate surface is however, the principal objective of the process. It is important that the deposition be uniform and without internal flaws.
ARTICLE IN PRESS 600
A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
The underlying flow and transport phenomena during CVD have a strong influence in the deposition layer thickness, its uniformity and purity. A mathematical model is essential for understanding the details of the CVD process and the influence of various operating parameters on the deposition rates. A variety of CVD geometries are presently available: the horizontal reactor, pan-cake reactor and barrel reactor being the most explored [1]. There is a continuing trend in microelectronics towards wafers of large diameter, higher quality, and a reduction in the minimum size of flaws. There is also a need to grow infrared sensors and radiation windows of various shapes and sizes. Thus, CVD technology calls for better understanding of the underlying transport phenomena and its relationship with the deposition characteristics. Numerical studies of the CVD process has been carried out by various authors to investigate the role of transport phenomena in the overall film growth characteristics of epitaxial silicon. Moffat and Jensen [2] studied the effect of three-dimensional flow phenomena on critical parameters, such as the deposition rate on the substrate surface. Klejin et al. [3] solved the flow, energy and species transport equations in two-dimensional form in an axisymmetric geometry, taking thermo-diffusion and multi-component diffusion into account. Oh et al. [4] developed a mathematical model for the epitaxial silicon growth in a pancake reactor where equations governing momentum, energy and species transport were numerically solved. Duverneuil and Couderc [5] presented a two-dimensional model for growth of silicon during low-pressure CVD. It involved solving the stream function–vorticity equations, with the representative inter-wafer region as the physical domain. Mahajan [6] has reviewed mathematical models for thin-film deposition by the CVD process. The modelling of low-pressure CVD reactors in the above investigations focused on the reactive substrates. The substrate is to be treated as chemically reactive, when growth of thin films in the micron-range is encountered. In the present work, the deposited material is of greater thick-
tube wall α=30 ο
peripheral jet (Zn+Ar) outflow plane deposition surface central jet (H2 S+Ar) y axis x
Fig. 1. Schematic drawing of a model CVD reactor with a convex substrate.
ness, being in the range of millimetres. In effect, the substrate can be modelled as chemically passive. The mass transfer boundary conditions for such a surface cannot be expressed in a conventional form by physical reasoning alone. In the present work, the uncertainty in the boundary condition has been dealt with by an approximation. A schematic view of the CVD reactor with a concave substrate is shown in Fig. 1 Convex and a concave-shaped substrate blocks facing the flow have been separately considered. The mathematical model developed in the present work assumes the flow and scalar fields to be axisymmetric.
2. Mathematical formulation The mixing and reaction of argon jets carrying hydrogen sulphide and zinc vapour, and the deposition of the product, namely zinc sulphide on the substrate surface have been numerically simulated in the present study. The mathematical model has been developed by employing the following assumptions: 1. The transport of the individual species is governed by the convection-diffusion equation with (a) the velocities independently computed, (b) the production rates determined by a singlestep chemical reaction, and (c) the reactants and products in argon constituting a dilute solution. The dilute solution approximation enables the higher-order effects of mass transfer to be neglected.
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
2. Zinc sulphide is treated as a gaseous phase, and is also subjected to the transport equation as applicable to the other species. 3. The fluid properties such as density and viscosity that determine the flow field are taken to be those of argon. 4. The reactor temperature during the deposition process is practically constant. Hence, the viscosity and density of the carrier gas are taken to be spatially constant. The assumption involving zinc sulphide requires an explanation. The product ZnS of the chemical reaction is composed of micron-sized particles. As the fluid forces exerted on these particles (proportional to the surface area) overwhelm their weight (proportional to the volume), they are transported in the carrier gas with the local gas velocity. Thus, zinc sulphide can be considered to be a passive scalar, driven by the base flow field towards the substrate, where it is deposited at a rate fixed by the mass-diffusion fluxes. Deposition occurs at the tube walls as well, but the focus of the present investigation is towards estimating deposition rates on the substrate held normal to the flow direction. The present study has been developed for lowpressure CVD reactors [7,8]. For a typical reactor, the following representative parameters can be assumed: diameter ¼ 0:1 m, bulk pressure ¼ 30 mbar, bulk temperature ¼ 900 K, characteristic flow rate of the inner jet ¼ 500 (standard) cm3 / min, the carrier gas being argon. The Reynolds number under these conditions is 100, the Prandtl number being of the order of unity. The thermophysical properties of the gases involved in the reaction are summarized in Table 1.
601
To assess the importance of buoyancy forces arising from spatial density differences, the Grashof number was calculated on the basis of the largest possible mass fractions of the respective species, the largest characteristic dimension, namely the radius, with other properties evaluated at the mean reactor pressure and temperature. The Grashof number was found to be 260, with the mixed convection parameter Gr=Re2 working out to be 0.026. This number is much smaller than unity, indicating that buoyancy can barely influence transport by forced convection. In addition to buoyancy, viscous dissipation is also negligible in the present conext. Thus, laminar incompressible flow conditions can be assumed in the analysis of the flow and scalar concentration fields. The thickness of the material deposited in the reactor is sufficiently small, being of the order of a few millimetres. Hence, the change in the shape and size of the flow space can be neglected during simulation. The continuity equation and the Navier–Stokes equation for laminar, incompressible flow are: r u ¼ 0;
ð1Þ
qu 1 þ r ðuuÞ ¼ rp þ r2 u: qt Re
ð2Þ
The governing equation for transport of the ith species involves diffusive mass fluxes Ji ; relative to mass-averaged velocity u; of the gas mixture. If oi is the mass fraction of the ith species, the associated transport equation is qðroi Þ þ r ðruoi Þ ¼ r Ji þ Si : qt
ð3Þ
Table 1 Transport properties of the constituent gases in a low-pressure CVD reactor at 917 K and 30 mbar
3
Density (kg/m ) Kinematic viscosity (m2 /s) Mass diffusivity in argon (m2 /s) Molecular mass (kg/kmol) kf (m3 /mol s)
Argon
H2 S
Zn vapour
H2
ZnS
0.0157 5.3E-03 — 39.94
0.0134 — 3.8E-03 34.08
0.0257 — 3.28E-03 65.38 106
7.9E-04 — 1.75E-2 2.016
0.0383 — 2.46E-3 97.38
Dashes indicate that these values are not required in the mathematical model.
ARTICLE IN PRESS 602
A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
Here the total diffusive mass flux Ji of species i is composed of Jci ; the flux due to Fickian diffusion caused by a concentration gradient, the pressure diffusion Jpi ; the forced diffusion Jfi ; and the thermal diffusion JTi : Hence Ji ¼ Jci þ Jpi þ Jfi þ JTi
ð4Þ
In a low-pressure CVD reactor the pressure gradients are small and thus pressure diffusion is not important. Gravitational effects are also not important. Without external forces, such as an electric field, forced diffusion is negligible. Of interest to the present work are CVD reactors that are temperature controlled in the sense that the temperature variation in the reactor is small. Hence, the thermo-diffusion effect can be neglected in the present formulation. Thus, in the dilute solution, only Fickian diffusion contributes to species diffusion. The possibility of thermo-diffusion effects in a low-pressure CVD reactor requires re-examination. The parameters employed in the numerical study were based on those employed in prototype reactors, for example the one reported by Zhenyi et al. [9]. The temperature levels in the reactor average to 900 K ð75 K), while the substrate temperature is slightly larger than 923 K, and is quite uniform. Thermo-diffusion effects refer to mass transfer set up by thermal gradients. It is a higher-order effect, in comparison to mass transfer by concentration gradients. Since the substrate and the reactor temperatures are quite close, it suggests, on an average that thermo-diffusion effects should be small. Also, the substrate is at a higher temperature than that of the reacting gases. Since heat release during the chemical reaction of interest is small, the product gas temperature is close to that of the reactants. The thermal gradients that are now set up between the substrate and the gases in the vicinity are small; further, they are of a sign that marginally reduces the deposition rate. In addition, the uniformity of temperature ensures that diffusion of zinc sulphide within the substrate is inhibited. Hence, the substrate may be assumed to be passive and thermo-diffusion effects to be negligible. The source term Si in Eq. (3) represents the creation or destruction of the ith species during
the chemical reaction. In the present study, the chemical reaction has been modelled as an elementary single-step, homogeneous reaction taking place in the gas phase. The emphasis of the present work is on deposits of the order of a few millimeters. The substrate surface reacts with the zinc sulphide particles initially when the deposited layer is a few microns thick. Hence, for later times, the substrate can be taken to be chemically inert. The gas-phase reaction is taken to be instantaneous, in the sense that the reaction constants are large. Since zinc sulphide, on formation, is a solid particle, the possibility of dissociation is small and the reaction has been considered to be irreversible. The reaction modelled in the present work contains four species and is expressed as follows: kf
nA A þ nB B 2 nC C þ nD D; kr
ð5Þ
where the n’s are the stoichiometric coefficients, and kf and kr are the forward and reverse reaction rates. The value of the forward reaction rate is given in Table 1. The net production rate RC for zinc sulphide (species ‘‘C’’ Eq. (5)) is given by RC ¼ RfC RrC ¼ kf ðroA ÞðroB Þ kr ðroC Þ ðroD Þ:
ð6Þ
The coefficient kr has been set to zero in the present work. The dimensionless form of the transport equation as applicable for the dilute solution of the species in argon is written as follows ([5,10]): qðroi Þ 1 Di;o þ r ðruoi Þ ¼ r ðrDi roi Þ qt Pem DA;o d r Mi Ri : ð7Þ þ ni oA;o c Uc Mc Mc Here, suffix ‘‘c’’ refers to the carrier gas. The notation used in Eq. (7) is defined in the nomenclature. 2.1. Boundary conditions Numerical calculations of the present work showed that the initial transients were short-lived, and quantities such as the deposition rates rapidly reached a steady state. Hence, results for the
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
steady fields alone have been discussed. These have been obtained by starting with zero initial conditions for the field variables over the entire reactor. Eqs. (1), (2) and (7) have been solved for the following boundary conditions: No-slip boundary conditions for velocity are prescribed at the walls and the substrate. The inlet velocities are fully prescribed at the jet inlets. The fully developed boundary condition corresponding to the zero normal derivative has been applied at the exit plane. The outflow plane condition is a source of uncertainty in the numerical calculations. To minimize the possibility of errors, particularly in the region between the inflow plane and the substrate, the outflow plane was located sufficently far away from the substrate. Numerical experiments showed that the deposition rates did not depend on the location of the outflow plane if the latter was more than a distance of one reactor diameter from the substrate. The species concentrations are fully prescribed at the jet locations. The zero normal derivative boundary conditions are prescribed at the solid walls and also at the substrate surface. As the species penetration at the walls is essentially by the diffusion fluxes, the zero normal derivative condition is equivalent to the walls are impermeable to the respective species. The outflow boundary condition is once again of the zero derivative boundary condition type. Previous researchers treated the substrates to be chemically active. Thus, the relevant boundary condition was of the mixed-type, incorporating a normal derivative term from the diffusion flux and a reaction term that is a function of the local concentration. For a passive surface this equation reduces to a zero normal derivative impermeability condition. A strictly impermeable surface would not permit any deposition of zinc sulphide at the substrate. It leads to the difficulty that in the presence of deposition, a correct boundary condition for ZnS is not available. To circumvent this difficulty, a special treatment was adopted for ZnS, based on the following reasoning: Deposition is driven by gradients in the mass fraction of ZnS, near the substrate. Under these conditions, the deposition rate for a clean surface is bound to be larger than that for a surface already carrying a
603
layer of deposit. For a clean and fresh substrate, the equivalent boundary condition for zinc sulphide is one of zero mass fraction. This condition has been enforced in the numerical computations. It is now understandable that the estimates from the simulation of the present study will provide an upper bound of the deposition rates. The Dirichlet boundary condition referred above is a special case of the Robin boundary condition applicable to a reactive surface. The generalized boundary condition in dimensionless form is
qoZnS ¼ BioZnS ; qn
where oZnS is the concentration of zinc sulphide, n is a unit outward drawn normal on the substrate, and Bi, the Biot number contains details of the surface reaction rates. The Dirichlet boundary condition is recovered in the limiting case when the dimensionless parameter Bi (and hence the rate constant) goes to infinity. The Neumann condition is realized for Bi ¼ 0; namely a non-reacting surface for which no deposition can take place. The boundary condition for deposition adopted in the present work is expected to work well for bulk growth of zinc sulphide in the polycrystalline form. It is free of any directional effects, and hence cannot impart a structure to the growing medium. In epitaxial growth, the layer is only a few microns thick and the passive substrate approximation is not expected to be valid. 2.2. Range of parameters Numerical calculations have been carried out for fluid flow and mass transfer over a range of parameters. The parameters include the Reynolds number based on the incoming velocity of the peripheral jets and the characteristic dimension, namely the tube diameter. For normal operating conditions, the typical value of Re was found to be 100 (Section 1). Results have been presented in the following sections for Reynolds numbers upto 200. Other dimensionless parameters are the block position Zb ; the ratio of the axial to the peripheral jet velocity Vr ; and the block size rb : The block thickness has not been neglected in the
ARTICLE IN PRESS 604
A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
calculations. In most runs, it was set to be equal to three cell thicknesses. The choice of three cells was motivated by the convenience in applying the no-slip condition. Since the internal grid was generated by solving an elliptic partial differential equation, the cell thicknesses for the concave and convex substartes are unequal. This factor does not affect the solutions upstream of the block. The dimensionless values employed in the simulation are: Vr ¼ 2; rb ¼ 0:40 and Zb ¼ 0:6–1. Reynolds number over the range 10–200 have been considered. The mass fractions of zinc and hydrogen sulphide on the inflow plane were selected to make the latter 20% in excess of the stoichiometric value. The binary diffusion coefficients with respect to argon were determined by the Chapman–Enskog formula derived from the kinetic theory of gases [11]. The Peclet numbers for the four species were found to be in the range 200– 400 corresponding to a Reynolds number of 100. The grid Courant number (rather than the grid Fourier number) was the controlling factor in fixing the permissible time step in the explicit time marching calculation. A reactor temperature of 900 K was chosen for simulation. This is quite similar to values of 900– 1000 K quoted by Zhenyi et al. [9] for a lowpressure CVD reactor for zinc sulphide. The change in the average temperature of the reactor changes the fluid properties, and most significantly the binary diffusion coefficients. These coefficients scale as T 1:5 : Hence, for an increase in temperature from 900 to 1000 K, the binary diffusion coefficient increases by a factor 1.17, and the Peclet number decreases by 15%. Since a wide Peclet number range has been covered, the property variation due to bulk temperature changes are covered by the parameter space of the study. The transport properties appearing in the mathematical formulation of Section 2 are summarized in Table 1.
solved using a finite volume formulation. Here, the governing equations are integrated over an arbitrary prism-shaped cell. The advantage of this approach is that with an appropriate grid generation algorithm, flow in complex geometries can be simulated. Specifically, there is no need for a coordinate transformation. Hence, flow in the tubular reactor can be computed without transforming the equations to a cylindrical coordinate system. A collocated grid in which all dependent variables are prescribed at identical nodal points was used. The numerical algorithm has been developed in three dimensions, and is applicable to internal as well as external flow geometries. The two-dimensional solution in x–y coordinates was derived from the original threedimensional code by considering only a sector of the tube carrying three cells in the angular direction. Two dimensionality could be enforced by applying the periodicity condition in velocity and concentration at the first and the third cells. For the axisymmetric geometry such as the reactor tube, changes in cell areas of the mesh with distance from the axis are fully accounted for during the grid generation procedure. The radial coordinate of the reactor shown in Fig. 1 is aligned with the y-coordinate of the computational domain. Body conforming grids have been generated for the present study by solving elliptic partial differential equations for the nodal coordinate functions xðx; ZÞ and yðx; ZÞ: These map the coordinates (x;y) in the physical domain to a uniform grid on a square in the (x;Z) plane. The flow and transport equations are solved by the finite volume method in the physical domain. The mapping can, however, be used to get uniform and orthogonal coordinate lines throughout the physical domain [12]. The grid generation equations utilized in the present work are q qx q 1 qx f þ ¼ 0; qx qx qZ f qZ
ð8Þ
q qy q 1 qy f þ ¼ 0; qx qx qZ f qZ
ð9Þ
3. Numerical procedure The governing differential equations of flow and heat transfer in Cartesian coordinates have been
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
where f is the distortion function given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2Z þ y2Z f ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2x þ y2x
ð10Þ
with the subscripts indicating partial derivatives. Eqs. (8) and (9) are solved with the boundary condition xx xZ þ yx yZ ¼ 0 which imposes orthogonality on the grid at the domain boundaries. Examples of grids generated in the present work are shown in Fig. 2. In the finite volume method, the physical domain is divided into a number of contiguous hexahedral control volumes whose vertices are generated by the grid generation technique. The variables are specified in the collocated manner at the centroids of the control volume. Momentum interpolation is used to avoid pressure–velocity decoupling in the numerical solution. The treatment adopted here follows closely that of Eswaran and Prakash [13] and Prabhakar et al. [14]. The species equations with the reaction terms are solved in a similar manner. An explicit time marching procedure, starting from quiescent initial conditions and specified boundary conditions is
605
used to obtain the steady state solution. With the dilute solution assumption, the species fields do not affect the flow field. Thus the Navier–Stokes and continuity equations are first solved till steady state, and the species equations are then solved using the steady velocity profiles. For the reaction of zinc with hydrogen sulphide, the forward reaction rate constant is quite large, and leads to numerical instabilities. To address this difficulty, an operator splitting approach is used for the species equations [15]. Here, the governing transport equation for each species is split into two parts and separately solved. In the first step one solves q 1 Di;o ðroi Þ þ r ðruoi Þ rðroi Þ qt Pem DA;o ¼ 0 over a time step Dt
ð11Þ to obtain the intermediate value oi ðt þ DtÞ from oi ðtÞ: The intermediate value is used as the initial condition to solve q ðroi Þ ¼ Ri over the time step Dt ð12Þ qt to obtain the final value of oi ðt þ DtÞ: This procedure is repeated for each time step. The first
r/D
0.5
0.25
0 0
0.5
1
1.5
2
1.5
2
x/D
(a)
r/D
0.5
0.25
0 0
(b)
0.5
1
x/D
Fig. 2. Finite volume grids used in the transport calculations for the reactor with (a) concave substrate and (b) convex substrate.
ARTICLE IN PRESS 606
A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
step is solved in the advection–diffusion finitevolume framework using an explicit scheme. The second step is solved by linearizing the reaction term and using an analytical solution. The advantage of operator splitting is that the high reaction rate constant does not impose an unduly small time-step restriction on the solver, as the reaction term is analytically treated. The present approach also led to an improvement in the overall computational speed. All calculations were carried out on a P-III, 600 MHz machine with 512 MB RAM. Typical CPU times required for completing a single run were found to be 10 h for the flow field and 6–8 days for the mass transfer calculations. Results obtained in the present study showed that the flow in the reactor is dominated by a recirculation pattern. Hence, there was a need to resolve not only the boundary layers, but also several regions away from the boundaries, the recirculation zone for example. A uniform grid was thus opted for, so that all regions in the reactor could be resolved without any bias. Based on a detailed grid independence study, a 83 48 grid in the x and y directions, respectively, was employed in all the computations. The time step was determined on the basis of CFL as well the Fourier number criteria. The computer code was validated against the experiments of Kaushik and Rubin [16], as well as the confined jet mixing solution presented by Seider and Churchill [17].
4. Results and discussion Results for fluid flow distribution and mass transfer rates in the CVD reactor are presented in this section. The flow field is shown by velocity vectors; relevant streamlines are superimposed to identify the recirculation pattern, specifically its size, strength and location. Contours of the steady state production rate of zinc sulphide have been shown to illustrate the size of the reaction zone. The mass fraction contours of zinc sulphide have also been shown with respect to the mass of the carrier gas entering the reactor. Subsequently, the deposition rates of zinc sulphide for the convex and the concave substrates have been compared.
4.1. Velocity field The flow patterns in CVD reactors with convex and concave substrates are shown in Figs. 3(a–f). While the Reynolds number is varied form 10 to 200, other quantities have been kept fixed in the simulation. Specifically, the velocity ratio of the central to the peripheral jet is 2, the injection angle of the peripheral jet is 30 towards the reactor axis, the blockage ratio of the substrate is 0.4 and the block distance from the jets is unity, with respect to the reactor diameter. Flow in the reactor with a concave substrate (left column of Fig. 3) contains a stagnation zone ahead of the substrate and a recirculation zone closer to the inflow plane. The latter is shown in Fig. 3 by solid streamline contours. The jets are seen to mix rapidly in terms of the distance along the axis. The reactor is filled with a characteristic flow pattern, that is quite independent of the inflow conditions. The recirculation pattern becomes stronger with an increase in Reynolds number. Simultaneously, the size of the stagnation zone ahead of the substrate decreases. The region adjacent to the reactor wall is also stagnant, and continues to be in this state for higher Reynolds numbers. The annular gap above the substrate through which the flow is exhausted is a region of high velocity. The right column of Fig. 3 shows the flow field in a reactor with a convex substrate for Reynolds numbers of 10, 100 and 200. A recirculation pattern above the peripheral jet is seen at Reynolds numbers of 100 and 200, but is generally weaker compared to a concave substrate. The region near the reactor wall has higher velocities. The stagnation zone ahead of the substrate is completely absent. These results are understandable, because the substrate is now partially streamlined with respect to the flow direction and offers less of an obstruction to the mean flow. 4.2. Species distribution Among the four species namely, zinc, hydrogen sulphide, zinc sulphide and hydrogen, the species of importance to the present study is zinc sulphide. The contours of the production rates of zinc sulphide have been determined from the fourth
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
607
Fig. 3. Velocity vectors for the reactor with a concave substrate: (a) Re ¼ 10; (b) Re ¼ 100; (c) Re ¼ 200; convex substrate: (d) Re ¼ 10; (e) Re ¼ 100; (f) Re ¼ 200: The solid line contour indicates the zone of reversed flow.
term of Eq. (7). They have been plotted as a function of Reynolds number in Fig. 4. The contours reflect the primary positions in the reactor where the chemical reaction is in progress. Since the size of the reaction zone alone is of interest, the dimensional production rates have been normalized between 0 and 10. Fig. 4 shows that the reaction is confined to the interfacial region between the central and the peripheral jets. The zone is aligned with the local flow direction. For a concave substrate (left column), the reaction zone exhibits a curvature oriented towards the annular gap. The curvature is smaller for the convex substrate (right column). The reaction zone is extended towards the substrate with an increasing on Reynolds number. These trends are clearly visible at a Reynolds number of 200. Fig. 4 shows
that the block position and shape affect the size of the reaction zone by influencing the flow field, but their direct influence on the chemical reactions is negligible. Specifically, the transport and deposition of zinc sulphide in the reactor follow the rules of a passive scalar. The mass fraction contours of zinc sulphide for the two substrates are presented in Fig. 5. The contour labels are in proportion to the magnitude of the mass fraction. For the concave substrate (left column of Fig. 5) the mass fractions close to it are small. The influence of increasing the Reynolds number is to stretch the contours in the flow direction. The consequent possibility of a higher deposition rate is, however, balanced by the prevalant stagnation region ahead of the substrate. For the convex substrate, an increase in Reynolds
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
608
0.25
0.5
1
r/D
r/D
0.5
0 0
0.5
(a)
r/D
r/D
1
x/D 0.5
1
0.25
0 0
0.5
1
0.25
0 0
1
x/D
(b)
0.5
1
x/D
(e)
0.5
0.5
0.25
10
5
r/D
r/D
0.5
(d)
0.5
(c)
3 1
0 0
1
x/D
0 0
0.25
1
0.5
0 0
1
x/D
0.25
(f )
7
1
5 3
0.5
1
x/D
Fig. 4. Contours of the production rates of zinc sulphide for the reactor with a concave substrate: (a) Re ¼ 10; (b) Re ¼ 100; (c) Re ¼ 200; convex substrate: (d) Re ¼ 10; (e) Re ¼ 100; (f) Re ¼ 200: Numbers indicate relative magnitudes.
number stretches the zinc sulphide contours towards the block (right column of Fig. 5). The absence of a stagnation region in this region suggests the possibility of higher deposition rates in this arrangement. For both geometries, Fig. 5 shows that a Reynolds number of 10 produces nearly vertical mass fraction contiurs. It is indicative of a stratification in the zinc sulphide concentration with low-concentration zones being formed in the vicinity of the substrates. From the viewpoint of uniformity of the thickness of the deposited layer, this is an advantage. The corresponding deposition rates computed as an average over the surface are expected to be low. 4.3. Wall shear stress and deposition rates It is of interest to see if the influence of the flow fields in Fig. 3 on mass transfer can be character-
ized in terms of the wall shear stresses. The shear stress distribution at the substrate surfaces are shown in Fig. 6 at a Reynolds number of 100, for two block positions. In view of the stagnation regions prevailing adjacent to the reactor tube as well as the substrate surface, these are regions of low-shear stress. At the substrate, the stress levels were seen in the computations to increase with increasing Reynolds number. At the reactor wall, the stresses were actually found to diminish with Reynolds number [18]. Fig. 6 shows that the shear stresses are higher for the convex substrate as compared to the concave. The tip of the substrate is a region of high stresses, where large velocities are encountered. A description in terms of the distribution of the wall shear stresses is useful for passive scalars such as thermal energy and non-reacting gases. Here the shear stresses can be empirically correlated by analogy to the wall fluxes of the scalar property.
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
0.25
0.5
7 9 5 10 10 3 9
5 3 5
3
1
r/D
r/D
0.5
7 9 9 5
0.25
7
0 0
1
5
3
0.5
3
5
7 1
r/D
1 r/D
0.5
1
5
3
5
9 9
7
0.25
5 1
1 0.5
(b)
0 0
1
0.5
1 x/D
0.5
3
1
0.25
9
9
5
7
1 r/D
r/D
0.5
(e)
x/D
7
9 5
0.25
3
0.5
9
0 0
1
x/D
(f )
9
5 7
3 9 1
1 0 0
1
3
3 0 0
1
x/D
9 7
3 1
7
5
0.5
7 9
0.25
1
3
(d)
x/D 0.5
0 0
1
(a)
(c)
609
3
5
7
0.5
1
x/D
Fig. 5. Contours of the mass fraction of zinc sulphide for the reactor with a concave substrate: (a) Re ¼ 10; (b) Re ¼ 100; (c) Re ¼ 200; convex substrate: (d) Re ¼ 10; (e) Re ¼ 100; (f) Re ¼ 200: Numbers indicate relative magnitudes.
The relationship between the fluxes and the stresses follows a direct proportionality rule, in the sense that regions of high shear are accompanied by high fluxes [19]. The correlation cannot be unconditionally extended to systems with reactions. The local wall fluxes of zinc sulphide are shown in Fig. 7. They are related to the material deposition rates at the curved blocks. The reactor wall and the substrate are regions of low deposition, while the edges of the substrate have high deposition rates of zinc sulphide. The wall shear stress and the wall fluxes of the transported species (and hence the deposition rates) thus show direct similarity. Hence, in the present work, the correlation between them is broadly realized under reaction conditions as well. The similarity between the wall shear stress and wall fluxes of zinc sulphide can be understood with reference to Fig. 4. Here it is revealed that the size
of the reaction zone is quite narrow and is located closer to the inflow plane of the jets. Hence, the deposition process in the reactor is dominated by scalar transport, for which the stress–flux similarity is applicable. A discussion on the relationship between shear stress and scalar wall fluxes for a reactor with a flat substrate has been presented by Saxena et al. [20]. Depositions of zinc sulphide would occur at all solid surfaces including the substrate and the tube walls. However, of interest here is the deposition at the substrate, and the effect of the surface shape on it. The local deposition rates for the convex and concave surfaces are presented in Fig. 7. It is clear that the deposition rates are higher for the convex substrate (Fig. 7b) when compared to the concave (Fig. 7a). The latter, however displays a greater uniformity in thickness, and only a marginal
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
610 0.4
0.4
0.35
0.35
concave
convex
0.3
radial coordinate
radial coordinate
0.3
0.25
0.2
0.15
0.25
0.2
0.15
0.1
0.1
Z b =0.6
0.05
Z b =0.6
0.05
Z b =1.0 o
Z b =1.0 o 0
0 0
0.05
0.1
(a)
0.15
0.2
0.25
0
0.3
0.2
0.4
0.6
(b)
Cf
0.8
1
1.2
1.4
Cf
Fig. 6. Wall shear stress profiles over the reactor substrates at Re ¼ 100; Vr ¼ 2; (a) concave, (b) convex.
0.4
0.4
0.35
0.35
0.25 0.2 0.15
Z b=1.0, Re=10 Z b=1.0, Re=100 Z b=1.0, Re=200 o Z b=0.6, Re=100
0.1 0.05 0
(a)
convex
0.3
radial coordinate
radial coordinate
concave
0
0.5
1
1.5
2
2.5
3
deposition rate (µm/min)
0.3 0.25 0.2 0.15
Z b=1.0, Re=10 Z b=1.0, Re=100 Z b=1.0, Re=200 o Z b=0.6, Re=100
0.1 0.05 0
3.5
(b)
0
0.5
1
1.5
2
2.5
3
deposition rate (µm/min)
Fig. 7. Deposition rate over the substrate as a function of Reynolds number; Vr ¼ 2; (a) concave, (b) convex.
influence of Reynolds number and the block position. For a convex substrate, the deposition rates decreases with an increase in Reynolds number. The deposition rate also decreases when the substrate is brought closer to the inflow plane (Fig. 7b, Zb =0.6). For Re ¼ 10; the deposition rate is reasonably high and the variation with the radial coordinate is uniform. The total material deposited can be calculated by averaging the deposition rate in the radial direction. The
deposited material is slightly lower at Re ¼ 10 compared to Re ¼ 100 and 200, though this is a minor disadvantage. For both geometries, the deposition rates near the edges of the substrate are high. This is because the bulk of zinc sulphide produced in the reactor bypasses the deposition surface through the annular gap. This excess material has an opportunity to stick to the corners of the substrate.
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
In an experiemental study, Zhenyi et al. [9] have discussed the performance of a CVD reactor for zinc sulphide. The operating conditions of the reactor in the study are quite similar to the present work. In the quoted study, the deposition rate over a collection of flat surfaces oriented parallel and perpendicular to the flow is 0.3–1 mm/day, namely 0.2–0.7 mm/min. In the present numerical study, the concave substrate yielded an average deposition rate of 0.5 mm/min at a Reynolds number of 100, while it was 1.027 mm/min for the convex substrate. These values are quite close to those quoted, despite differences in the substrate geometry and orientation. The average deposition rate at a prescribed H2 S loading was computed for a convex substrate to be 1.027 mm/min compared to 0.542 mm/min for a concave substrate, an improvement of a factor of two, at a Reynolds number of 100. The doubling of the deposition rate is a significant result. The increase in the deposition rate for a convex substrate with respect to the concave is due to the following factors: (i) a better streamlining of the flow field (Fig. 3), (ii) a reduction in the size of the stagnation zone (Fig. 3), and a consequent increase in the size of the reaction zone (Figs. 4 and 5). For the concave substrate, the size of the stagnation zone barely changes with Reynolds number. Hence, the deposition rates (as seen in Fig. 7) are lower and more uniform as compared to the convex substrate, and are insensitive to Reynolds number.
5. Conclusions Fluid flow patterns and mass transfer phenomena inside a low-pressure CVD reactor have been studied via a mathematical model. Zinc vapour and hydrogen sulphide, injected along with argon react chemically in the reactor to form zinc sulphide. Of special interest is the deposition characteristics of zinc sulphide on a passive substrate held normal to the direction of the incoming flow. The following conclusions have been arrived at in the present work: 1. When a concave substrate is employed, a large stagnation region is formed where the fluid
611
velocities are small. A considerable fraction of the product species bypasses the substrate. Hence-the deposition rate is small but uniform. The deposition rate increases marginally with Reynolds number because of a slight reduction in the size of the stagnation zone. 2. A convex substrate is better streamlined, and allows the bulk of flow and the reacting species to come closer to it. This results in a smaller stagnation zone, whose size depends on Reynolds number. Species advection by the carrier gas, in addition to species diffusion significantly contributes towards the deposition of zinc sulphide on the substrate. Hence, the deposition rate increases over and above what is achieved by the concave substrate. 3. Both convex and concave substrates share unfavourable edge effects. 4. Positioning the substrate closer to the inflow plane increases the size of the stagnation zone and hence diminishes the deposition rates for both geometries.
References [1] M.L. Hitchman, K.F. Jensen, Chemical vapour Deposition, Academic Press, New york, 1993. [2] H.K. Moffat, K.F. Jensen, J. Electrochem. Soc. 135 (2) (1988) 459. [3] C.R. Klejin, Th. H. van der Meer, C.J. Hoogendoorn, J. Electrochem. Soc. 136 (11) (1989) 3423. [4] I.H. Oh, C.G. Takoudis, G.W. Neudeck, J. Electrochem. Soc. 138 (2) (1991) 554. [5] P. Duverneuil, J.P. Couderc, J. Electrochem. Soc. 139 (1) (1992) 296. [6] R.L. Mahajan, Adv. Heat Transfer 28 (1996) 339. [7] C.E. Morosanu, Thin Films by Chemical vapour Deposition, Elsevier, Amsterdam, 1990. [8] S.A. Campbell, The Science and Engineering of Microelectronic Fabrication, Oxford University Press, Oxford, UK, 1996. [9] F. Zhenyi, C. Yichao, H. Yongliang, Y. Yaoyuan, D. Yanping, Y. Zewu, T. Hongchang, X. Hongtao, W. Heming, J. Crystal Growth 237–239 (2002) 1707. [10] O. Levenspiel, Chemical Reaction Engineering, Wiley, Newyork, 1999. [11] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, Wiley, New York, 1960. [12] J.F. Thompson, Z.U.A. Warsi, C. Wayne Mastin, Numerical Grid Generation, North-Holland, Amsterdam, 1985.
ARTICLE IN PRESS 612
A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
[13] V. Eswaran, S. Prakash, A finite volume method for Navier–Stokes equations, Proceedings of the Third Asian CFD Conference, Vol. 1, Bangalore, 1998. [14] V. Prabhakar, G. Biswas, V. Eswaran, Numer. Heat Transfer 41 (Part A) 41 (2002) 1. [15] K. Muralidhar, M. Verghese, K.M. Pillai, Numer. Heat Transfer, 23 (Part A) (1993) 99. [16] S. Kaushik, S.G. Rubin, Comput. Fluids 24 (1) (1995) 27.
[17] W.D. Seider, S.W. Churchill, AIChE J. 17 (3) (1971) 704. [18] A.K. De, Modeling of transport phenomena in a low pressure CVD reactor, Master’s Thesis, IIT Kanpur, India, 2002. [19] W.M. Kays, M.E. Crawford, Convective Heat and Mass Transfer, International Edition, McGraw-Hill, New York, 1993. [20] V. Saxena, K. Muralidhar, V. Eswaran, Sadhana 27 (6) (2002) 1.
Journal of Crystal Growth 267 (2004) 598–612
Modelling of transport phenomena in a low-pressure CVD reactor Arnab K. Dea, K. Muralidhara,*, V. Eswarana, V.K Wadhawanb a
Department of Mechanical Engineering, Indian Institute of Technology, IIT Kanpur, Kanpur (UP) 208 016, India b Laser Materials Division, Center for Advanced Technology, Indore 452 013, India Received 13 April 2004; accepted 16 April 2004
Communicated by J.J. Derby
Abstract The present work reports the analysis of transport phenomena in a low-pressure chemical vapour deposition (CVD) reactor. A two-dimensional unsteady model has been used to compute flow, chemical reactions and mass transfer rates within the reactor. The formation of zinc sulphide from the reaction of zinc vapour and gaseous hydrogen sulphide has been analysed. Two geometric designs of the reactor have been studied for assessing the overall deposition characteristics. The mathematical model comprises of the conservation equations of mass, momentum and the mass fractions of the reactant and product species. The governing partial differential equations have been solved by the finite volume method. The deposition surface has been taken to be chemically passive. There is an uncertainty in the boundary condition for the depositing species on the passive surface. The difficulty has been resolved by utilizing a homogenous Dirichlet condition at the surface. The present approach yields an upper bound on the deposition rate. The numerical results obtained show that the species concentration field, correlates closely to the flow field in the reactor. Over the range of Reynolds numbers studied namely 10–200, a recirculation pattern adjacent to the side walls is revealed and a stagnation zone forms ahead of the block. Reactant species mix predominantly in the zone close to the inflow plane and the product species are transported by advection. The reactor with a convex substrate streamlines the flow better as compared to the concave. As a result, higher concentrations of the product species reach the substrate. At a Reynolds number of 100, the deposition rate for the convex substrate was seen to increase by a factor of two with respect to the concave. r 2004 Elsevier B.V. All rights reserved. PACS: 81.10.Bk; 81.15.Gh; 02.60.Cb; 02.70.Fj; 51.10.ty; 51.20.td Keywords: A1. Computer simulation; A1. Convection; A1. Diffusion; A1. Fluid flows; A1. Mass transfer; A2. Growth from vapor
*Corresponding author. Tel.: +91-5122597182; fax: +91-5122597408. E-mail addresses: [email protected] (A. K. De), [email protected] (K. Muralidhar), [email protected] (V. Eswaran), [email protected] (V.K. Wadhawan). 0022-0248/$ - see front matter r 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2004.04.036
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
599
Nomenclature Cf Di D Mi Mc Pem Re rb kf kr s Uc Vr x; y; z
Zb
skin friction coefficient, tw =12rc Uc2 binary diffusion coefficient of species i in argon, m2 =s diameter of the reactor molecular weight of the ith species molecular weight of the carrier gas, namely argon Peclet number associated with mass transfer, Uc D=Di reynolds number, rUc D=m dimensionless radius of the substrate forward reaction rate, m3 /mol s backward reaction rate, m3 /mol s thickness of the deposited zinc sulphide, mm velocity of the central jet, m/s central-to-peripheral jet velocity ratio Cartesian coordinates in the physical domain with x pointing in the flow direction dimensionless distance of the block from the inflow plane
1. Introduction Chemical vapour deposition (CVD) is an industrial process for growing thin layers (0.01– 10 mm) of semiconductor materials, such as epitaxial silicon and gallium arsenide. Thin films are technologically important, especially in the field of microelectronics and as protective coatings for optical elements. As against thin films, the application of CVD to the bulk growth of monolithic compounds such as zinc sulphide and cadmium telluride has been considered in the present work. Infrared sensors are an important application where such materials are required. The deposited material has thicknesses in the range 1–10 mm. The reactor configuration specifically analysed operates in the batch-mode, and employs low-gas pressures and high temperatures. The CVD reactor is a large tube with inflow on one side and a deposition surface on the other. It contains nozzles at the inlet plane for separate
Greek letters a injection angle of the peripheral jet m dynamic viscosity of argon, Pa s stoichiometric coefficient for ith species ni in the chemical reaction oi mass fraction of ith species oA;o mass fraction of hydrogen sulphide at the inlet Rfj forward rate of reaction of jth species, mol m3 s r Rj backward rate of reaction of jth species, mol m3 s rc density of the carrier gas, namely argon, kg/m3 rZnS density of zinc sulphide, kg/m3 tw wall shear stress, N/m2 Subscript o inflow plane Superscripts g gas phase s surface properties
injection of the two reactants, each mixed with an inert carrier gas. The substrate is a flat or a curved surface that is placed at normal incidence to the incoming flow. The flow bypass occurs at the periphery of the substrate. The incoming gases chemically react, resulting in fine particles of the product. These are carried by the main flow and deposited on the substrate surface. The particles, in addition may be adsorbed and later mobilized by diffusion within the substrate. The other gaseous products are exhausted along with the main flow from the reactor. In the present study, the reactants are zinc vapour and hydrogen sulphide, while the carrier gas is argon. The product species of importance to the present work is zinc sulphide. Deposition of zinc sulphide can occur at all the solid surfaces of the reactor. The deposition on the substrate surface is however, the principal objective of the process. It is important that the deposition be uniform and without internal flaws.
ARTICLE IN PRESS 600
A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
The underlying flow and transport phenomena during CVD have a strong influence in the deposition layer thickness, its uniformity and purity. A mathematical model is essential for understanding the details of the CVD process and the influence of various operating parameters on the deposition rates. A variety of CVD geometries are presently available: the horizontal reactor, pan-cake reactor and barrel reactor being the most explored [1]. There is a continuing trend in microelectronics towards wafers of large diameter, higher quality, and a reduction in the minimum size of flaws. There is also a need to grow infrared sensors and radiation windows of various shapes and sizes. Thus, CVD technology calls for better understanding of the underlying transport phenomena and its relationship with the deposition characteristics. Numerical studies of the CVD process has been carried out by various authors to investigate the role of transport phenomena in the overall film growth characteristics of epitaxial silicon. Moffat and Jensen [2] studied the effect of three-dimensional flow phenomena on critical parameters, such as the deposition rate on the substrate surface. Klejin et al. [3] solved the flow, energy and species transport equations in two-dimensional form in an axisymmetric geometry, taking thermo-diffusion and multi-component diffusion into account. Oh et al. [4] developed a mathematical model for the epitaxial silicon growth in a pancake reactor where equations governing momentum, energy and species transport were numerically solved. Duverneuil and Couderc [5] presented a two-dimensional model for growth of silicon during low-pressure CVD. It involved solving the stream function–vorticity equations, with the representative inter-wafer region as the physical domain. Mahajan [6] has reviewed mathematical models for thin-film deposition by the CVD process. The modelling of low-pressure CVD reactors in the above investigations focused on the reactive substrates. The substrate is to be treated as chemically reactive, when growth of thin films in the micron-range is encountered. In the present work, the deposited material is of greater thick-
tube wall α=30 ο
peripheral jet (Zn+Ar) outflow plane deposition surface central jet (H2 S+Ar) y axis x
Fig. 1. Schematic drawing of a model CVD reactor with a convex substrate.
ness, being in the range of millimetres. In effect, the substrate can be modelled as chemically passive. The mass transfer boundary conditions for such a surface cannot be expressed in a conventional form by physical reasoning alone. In the present work, the uncertainty in the boundary condition has been dealt with by an approximation. A schematic view of the CVD reactor with a concave substrate is shown in Fig. 1 Convex and a concave-shaped substrate blocks facing the flow have been separately considered. The mathematical model developed in the present work assumes the flow and scalar fields to be axisymmetric.
2. Mathematical formulation The mixing and reaction of argon jets carrying hydrogen sulphide and zinc vapour, and the deposition of the product, namely zinc sulphide on the substrate surface have been numerically simulated in the present study. The mathematical model has been developed by employing the following assumptions: 1. The transport of the individual species is governed by the convection-diffusion equation with (a) the velocities independently computed, (b) the production rates determined by a singlestep chemical reaction, and (c) the reactants and products in argon constituting a dilute solution. The dilute solution approximation enables the higher-order effects of mass transfer to be neglected.
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
2. Zinc sulphide is treated as a gaseous phase, and is also subjected to the transport equation as applicable to the other species. 3. The fluid properties such as density and viscosity that determine the flow field are taken to be those of argon. 4. The reactor temperature during the deposition process is practically constant. Hence, the viscosity and density of the carrier gas are taken to be spatially constant. The assumption involving zinc sulphide requires an explanation. The product ZnS of the chemical reaction is composed of micron-sized particles. As the fluid forces exerted on these particles (proportional to the surface area) overwhelm their weight (proportional to the volume), they are transported in the carrier gas with the local gas velocity. Thus, zinc sulphide can be considered to be a passive scalar, driven by the base flow field towards the substrate, where it is deposited at a rate fixed by the mass-diffusion fluxes. Deposition occurs at the tube walls as well, but the focus of the present investigation is towards estimating deposition rates on the substrate held normal to the flow direction. The present study has been developed for lowpressure CVD reactors [7,8]. For a typical reactor, the following representative parameters can be assumed: diameter ¼ 0:1 m, bulk pressure ¼ 30 mbar, bulk temperature ¼ 900 K, characteristic flow rate of the inner jet ¼ 500 (standard) cm3 / min, the carrier gas being argon. The Reynolds number under these conditions is 100, the Prandtl number being of the order of unity. The thermophysical properties of the gases involved in the reaction are summarized in Table 1.
601
To assess the importance of buoyancy forces arising from spatial density differences, the Grashof number was calculated on the basis of the largest possible mass fractions of the respective species, the largest characteristic dimension, namely the radius, with other properties evaluated at the mean reactor pressure and temperature. The Grashof number was found to be 260, with the mixed convection parameter Gr=Re2 working out to be 0.026. This number is much smaller than unity, indicating that buoyancy can barely influence transport by forced convection. In addition to buoyancy, viscous dissipation is also negligible in the present conext. Thus, laminar incompressible flow conditions can be assumed in the analysis of the flow and scalar concentration fields. The thickness of the material deposited in the reactor is sufficiently small, being of the order of a few millimetres. Hence, the change in the shape and size of the flow space can be neglected during simulation. The continuity equation and the Navier–Stokes equation for laminar, incompressible flow are: r u ¼ 0;
ð1Þ
qu 1 þ r ðuuÞ ¼ rp þ r2 u: qt Re
ð2Þ
The governing equation for transport of the ith species involves diffusive mass fluxes Ji ; relative to mass-averaged velocity u; of the gas mixture. If oi is the mass fraction of the ith species, the associated transport equation is qðroi Þ þ r ðruoi Þ ¼ r Ji þ Si : qt
ð3Þ
Table 1 Transport properties of the constituent gases in a low-pressure CVD reactor at 917 K and 30 mbar
3
Density (kg/m ) Kinematic viscosity (m2 /s) Mass diffusivity in argon (m2 /s) Molecular mass (kg/kmol) kf (m3 /mol s)
Argon
H2 S
Zn vapour
H2
ZnS
0.0157 5.3E-03 — 39.94
0.0134 — 3.8E-03 34.08
0.0257 — 3.28E-03 65.38 106
7.9E-04 — 1.75E-2 2.016
0.0383 — 2.46E-3 97.38
Dashes indicate that these values are not required in the mathematical model.
ARTICLE IN PRESS 602
A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
Here the total diffusive mass flux Ji of species i is composed of Jci ; the flux due to Fickian diffusion caused by a concentration gradient, the pressure diffusion Jpi ; the forced diffusion Jfi ; and the thermal diffusion JTi : Hence Ji ¼ Jci þ Jpi þ Jfi þ JTi
ð4Þ
In a low-pressure CVD reactor the pressure gradients are small and thus pressure diffusion is not important. Gravitational effects are also not important. Without external forces, such as an electric field, forced diffusion is negligible. Of interest to the present work are CVD reactors that are temperature controlled in the sense that the temperature variation in the reactor is small. Hence, the thermo-diffusion effect can be neglected in the present formulation. Thus, in the dilute solution, only Fickian diffusion contributes to species diffusion. The possibility of thermo-diffusion effects in a low-pressure CVD reactor requires re-examination. The parameters employed in the numerical study were based on those employed in prototype reactors, for example the one reported by Zhenyi et al. [9]. The temperature levels in the reactor average to 900 K ð75 K), while the substrate temperature is slightly larger than 923 K, and is quite uniform. Thermo-diffusion effects refer to mass transfer set up by thermal gradients. It is a higher-order effect, in comparison to mass transfer by concentration gradients. Since the substrate and the reactor temperatures are quite close, it suggests, on an average that thermo-diffusion effects should be small. Also, the substrate is at a higher temperature than that of the reacting gases. Since heat release during the chemical reaction of interest is small, the product gas temperature is close to that of the reactants. The thermal gradients that are now set up between the substrate and the gases in the vicinity are small; further, they are of a sign that marginally reduces the deposition rate. In addition, the uniformity of temperature ensures that diffusion of zinc sulphide within the substrate is inhibited. Hence, the substrate may be assumed to be passive and thermo-diffusion effects to be negligible. The source term Si in Eq. (3) represents the creation or destruction of the ith species during
the chemical reaction. In the present study, the chemical reaction has been modelled as an elementary single-step, homogeneous reaction taking place in the gas phase. The emphasis of the present work is on deposits of the order of a few millimeters. The substrate surface reacts with the zinc sulphide particles initially when the deposited layer is a few microns thick. Hence, for later times, the substrate can be taken to be chemically inert. The gas-phase reaction is taken to be instantaneous, in the sense that the reaction constants are large. Since zinc sulphide, on formation, is a solid particle, the possibility of dissociation is small and the reaction has been considered to be irreversible. The reaction modelled in the present work contains four species and is expressed as follows: kf
nA A þ nB B 2 nC C þ nD D; kr
ð5Þ
where the n’s are the stoichiometric coefficients, and kf and kr are the forward and reverse reaction rates. The value of the forward reaction rate is given in Table 1. The net production rate RC for zinc sulphide (species ‘‘C’’ Eq. (5)) is given by RC ¼ RfC RrC ¼ kf ðroA ÞðroB Þ kr ðroC Þ ðroD Þ:
ð6Þ
The coefficient kr has been set to zero in the present work. The dimensionless form of the transport equation as applicable for the dilute solution of the species in argon is written as follows ([5,10]): qðroi Þ 1 Di;o þ r ðruoi Þ ¼ r ðrDi roi Þ qt Pem DA;o d r Mi Ri : ð7Þ þ ni oA;o c Uc Mc Mc Here, suffix ‘‘c’’ refers to the carrier gas. The notation used in Eq. (7) is defined in the nomenclature. 2.1. Boundary conditions Numerical calculations of the present work showed that the initial transients were short-lived, and quantities such as the deposition rates rapidly reached a steady state. Hence, results for the
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
steady fields alone have been discussed. These have been obtained by starting with zero initial conditions for the field variables over the entire reactor. Eqs. (1), (2) and (7) have been solved for the following boundary conditions: No-slip boundary conditions for velocity are prescribed at the walls and the substrate. The inlet velocities are fully prescribed at the jet inlets. The fully developed boundary condition corresponding to the zero normal derivative has been applied at the exit plane. The outflow plane condition is a source of uncertainty in the numerical calculations. To minimize the possibility of errors, particularly in the region between the inflow plane and the substrate, the outflow plane was located sufficently far away from the substrate. Numerical experiments showed that the deposition rates did not depend on the location of the outflow plane if the latter was more than a distance of one reactor diameter from the substrate. The species concentrations are fully prescribed at the jet locations. The zero normal derivative boundary conditions are prescribed at the solid walls and also at the substrate surface. As the species penetration at the walls is essentially by the diffusion fluxes, the zero normal derivative condition is equivalent to the walls are impermeable to the respective species. The outflow boundary condition is once again of the zero derivative boundary condition type. Previous researchers treated the substrates to be chemically active. Thus, the relevant boundary condition was of the mixed-type, incorporating a normal derivative term from the diffusion flux and a reaction term that is a function of the local concentration. For a passive surface this equation reduces to a zero normal derivative impermeability condition. A strictly impermeable surface would not permit any deposition of zinc sulphide at the substrate. It leads to the difficulty that in the presence of deposition, a correct boundary condition for ZnS is not available. To circumvent this difficulty, a special treatment was adopted for ZnS, based on the following reasoning: Deposition is driven by gradients in the mass fraction of ZnS, near the substrate. Under these conditions, the deposition rate for a clean surface is bound to be larger than that for a surface already carrying a
603
layer of deposit. For a clean and fresh substrate, the equivalent boundary condition for zinc sulphide is one of zero mass fraction. This condition has been enforced in the numerical computations. It is now understandable that the estimates from the simulation of the present study will provide an upper bound of the deposition rates. The Dirichlet boundary condition referred above is a special case of the Robin boundary condition applicable to a reactive surface. The generalized boundary condition in dimensionless form is
qoZnS ¼ BioZnS ; qn
where oZnS is the concentration of zinc sulphide, n is a unit outward drawn normal on the substrate, and Bi, the Biot number contains details of the surface reaction rates. The Dirichlet boundary condition is recovered in the limiting case when the dimensionless parameter Bi (and hence the rate constant) goes to infinity. The Neumann condition is realized for Bi ¼ 0; namely a non-reacting surface for which no deposition can take place. The boundary condition for deposition adopted in the present work is expected to work well for bulk growth of zinc sulphide in the polycrystalline form. It is free of any directional effects, and hence cannot impart a structure to the growing medium. In epitaxial growth, the layer is only a few microns thick and the passive substrate approximation is not expected to be valid. 2.2. Range of parameters Numerical calculations have been carried out for fluid flow and mass transfer over a range of parameters. The parameters include the Reynolds number based on the incoming velocity of the peripheral jets and the characteristic dimension, namely the tube diameter. For normal operating conditions, the typical value of Re was found to be 100 (Section 1). Results have been presented in the following sections for Reynolds numbers upto 200. Other dimensionless parameters are the block position Zb ; the ratio of the axial to the peripheral jet velocity Vr ; and the block size rb : The block thickness has not been neglected in the
ARTICLE IN PRESS 604
A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
calculations. In most runs, it was set to be equal to three cell thicknesses. The choice of three cells was motivated by the convenience in applying the no-slip condition. Since the internal grid was generated by solving an elliptic partial differential equation, the cell thicknesses for the concave and convex substartes are unequal. This factor does not affect the solutions upstream of the block. The dimensionless values employed in the simulation are: Vr ¼ 2; rb ¼ 0:40 and Zb ¼ 0:6–1. Reynolds number over the range 10–200 have been considered. The mass fractions of zinc and hydrogen sulphide on the inflow plane were selected to make the latter 20% in excess of the stoichiometric value. The binary diffusion coefficients with respect to argon were determined by the Chapman–Enskog formula derived from the kinetic theory of gases [11]. The Peclet numbers for the four species were found to be in the range 200– 400 corresponding to a Reynolds number of 100. The grid Courant number (rather than the grid Fourier number) was the controlling factor in fixing the permissible time step in the explicit time marching calculation. A reactor temperature of 900 K was chosen for simulation. This is quite similar to values of 900– 1000 K quoted by Zhenyi et al. [9] for a lowpressure CVD reactor for zinc sulphide. The change in the average temperature of the reactor changes the fluid properties, and most significantly the binary diffusion coefficients. These coefficients scale as T 1:5 : Hence, for an increase in temperature from 900 to 1000 K, the binary diffusion coefficient increases by a factor 1.17, and the Peclet number decreases by 15%. Since a wide Peclet number range has been covered, the property variation due to bulk temperature changes are covered by the parameter space of the study. The transport properties appearing in the mathematical formulation of Section 2 are summarized in Table 1.
solved using a finite volume formulation. Here, the governing equations are integrated over an arbitrary prism-shaped cell. The advantage of this approach is that with an appropriate grid generation algorithm, flow in complex geometries can be simulated. Specifically, there is no need for a coordinate transformation. Hence, flow in the tubular reactor can be computed without transforming the equations to a cylindrical coordinate system. A collocated grid in which all dependent variables are prescribed at identical nodal points was used. The numerical algorithm has been developed in three dimensions, and is applicable to internal as well as external flow geometries. The two-dimensional solution in x–y coordinates was derived from the original threedimensional code by considering only a sector of the tube carrying three cells in the angular direction. Two dimensionality could be enforced by applying the periodicity condition in velocity and concentration at the first and the third cells. For the axisymmetric geometry such as the reactor tube, changes in cell areas of the mesh with distance from the axis are fully accounted for during the grid generation procedure. The radial coordinate of the reactor shown in Fig. 1 is aligned with the y-coordinate of the computational domain. Body conforming grids have been generated for the present study by solving elliptic partial differential equations for the nodal coordinate functions xðx; ZÞ and yðx; ZÞ: These map the coordinates (x;y) in the physical domain to a uniform grid on a square in the (x;Z) plane. The flow and transport equations are solved by the finite volume method in the physical domain. The mapping can, however, be used to get uniform and orthogonal coordinate lines throughout the physical domain [12]. The grid generation equations utilized in the present work are q qx q 1 qx f þ ¼ 0; qx qx qZ f qZ
ð8Þ
q qy q 1 qy f þ ¼ 0; qx qx qZ f qZ
ð9Þ
3. Numerical procedure The governing differential equations of flow and heat transfer in Cartesian coordinates have been
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
where f is the distortion function given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2Z þ y2Z f ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2x þ y2x
ð10Þ
with the subscripts indicating partial derivatives. Eqs. (8) and (9) are solved with the boundary condition xx xZ þ yx yZ ¼ 0 which imposes orthogonality on the grid at the domain boundaries. Examples of grids generated in the present work are shown in Fig. 2. In the finite volume method, the physical domain is divided into a number of contiguous hexahedral control volumes whose vertices are generated by the grid generation technique. The variables are specified in the collocated manner at the centroids of the control volume. Momentum interpolation is used to avoid pressure–velocity decoupling in the numerical solution. The treatment adopted here follows closely that of Eswaran and Prakash [13] and Prabhakar et al. [14]. The species equations with the reaction terms are solved in a similar manner. An explicit time marching procedure, starting from quiescent initial conditions and specified boundary conditions is
605
used to obtain the steady state solution. With the dilute solution assumption, the species fields do not affect the flow field. Thus the Navier–Stokes and continuity equations are first solved till steady state, and the species equations are then solved using the steady velocity profiles. For the reaction of zinc with hydrogen sulphide, the forward reaction rate constant is quite large, and leads to numerical instabilities. To address this difficulty, an operator splitting approach is used for the species equations [15]. Here, the governing transport equation for each species is split into two parts and separately solved. In the first step one solves q 1 Di;o ðroi Þ þ r ðruoi Þ rðroi Þ qt Pem DA;o ¼ 0 over a time step Dt
ð11Þ to obtain the intermediate value oi ðt þ DtÞ from oi ðtÞ: The intermediate value is used as the initial condition to solve q ðroi Þ ¼ Ri over the time step Dt ð12Þ qt to obtain the final value of oi ðt þ DtÞ: This procedure is repeated for each time step. The first
r/D
0.5
0.25
0 0
0.5
1
1.5
2
1.5
2
x/D
(a)
r/D
0.5
0.25
0 0
(b)
0.5
1
x/D
Fig. 2. Finite volume grids used in the transport calculations for the reactor with (a) concave substrate and (b) convex substrate.
ARTICLE IN PRESS 606
A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
step is solved in the advection–diffusion finitevolume framework using an explicit scheme. The second step is solved by linearizing the reaction term and using an analytical solution. The advantage of operator splitting is that the high reaction rate constant does not impose an unduly small time-step restriction on the solver, as the reaction term is analytically treated. The present approach also led to an improvement in the overall computational speed. All calculations were carried out on a P-III, 600 MHz machine with 512 MB RAM. Typical CPU times required for completing a single run were found to be 10 h for the flow field and 6–8 days for the mass transfer calculations. Results obtained in the present study showed that the flow in the reactor is dominated by a recirculation pattern. Hence, there was a need to resolve not only the boundary layers, but also several regions away from the boundaries, the recirculation zone for example. A uniform grid was thus opted for, so that all regions in the reactor could be resolved without any bias. Based on a detailed grid independence study, a 83 48 grid in the x and y directions, respectively, was employed in all the computations. The time step was determined on the basis of CFL as well the Fourier number criteria. The computer code was validated against the experiments of Kaushik and Rubin [16], as well as the confined jet mixing solution presented by Seider and Churchill [17].
4. Results and discussion Results for fluid flow distribution and mass transfer rates in the CVD reactor are presented in this section. The flow field is shown by velocity vectors; relevant streamlines are superimposed to identify the recirculation pattern, specifically its size, strength and location. Contours of the steady state production rate of zinc sulphide have been shown to illustrate the size of the reaction zone. The mass fraction contours of zinc sulphide have also been shown with respect to the mass of the carrier gas entering the reactor. Subsequently, the deposition rates of zinc sulphide for the convex and the concave substrates have been compared.
4.1. Velocity field The flow patterns in CVD reactors with convex and concave substrates are shown in Figs. 3(a–f). While the Reynolds number is varied form 10 to 200, other quantities have been kept fixed in the simulation. Specifically, the velocity ratio of the central to the peripheral jet is 2, the injection angle of the peripheral jet is 30 towards the reactor axis, the blockage ratio of the substrate is 0.4 and the block distance from the jets is unity, with respect to the reactor diameter. Flow in the reactor with a concave substrate (left column of Fig. 3) contains a stagnation zone ahead of the substrate and a recirculation zone closer to the inflow plane. The latter is shown in Fig. 3 by solid streamline contours. The jets are seen to mix rapidly in terms of the distance along the axis. The reactor is filled with a characteristic flow pattern, that is quite independent of the inflow conditions. The recirculation pattern becomes stronger with an increase in Reynolds number. Simultaneously, the size of the stagnation zone ahead of the substrate decreases. The region adjacent to the reactor wall is also stagnant, and continues to be in this state for higher Reynolds numbers. The annular gap above the substrate through which the flow is exhausted is a region of high velocity. The right column of Fig. 3 shows the flow field in a reactor with a convex substrate for Reynolds numbers of 10, 100 and 200. A recirculation pattern above the peripheral jet is seen at Reynolds numbers of 100 and 200, but is generally weaker compared to a concave substrate. The region near the reactor wall has higher velocities. The stagnation zone ahead of the substrate is completely absent. These results are understandable, because the substrate is now partially streamlined with respect to the flow direction and offers less of an obstruction to the mean flow. 4.2. Species distribution Among the four species namely, zinc, hydrogen sulphide, zinc sulphide and hydrogen, the species of importance to the present study is zinc sulphide. The contours of the production rates of zinc sulphide have been determined from the fourth
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
607
Fig. 3. Velocity vectors for the reactor with a concave substrate: (a) Re ¼ 10; (b) Re ¼ 100; (c) Re ¼ 200; convex substrate: (d) Re ¼ 10; (e) Re ¼ 100; (f) Re ¼ 200: The solid line contour indicates the zone of reversed flow.
term of Eq. (7). They have been plotted as a function of Reynolds number in Fig. 4. The contours reflect the primary positions in the reactor where the chemical reaction is in progress. Since the size of the reaction zone alone is of interest, the dimensional production rates have been normalized between 0 and 10. Fig. 4 shows that the reaction is confined to the interfacial region between the central and the peripheral jets. The zone is aligned with the local flow direction. For a concave substrate (left column), the reaction zone exhibits a curvature oriented towards the annular gap. The curvature is smaller for the convex substrate (right column). The reaction zone is extended towards the substrate with an increasing on Reynolds number. These trends are clearly visible at a Reynolds number of 200. Fig. 4 shows
that the block position and shape affect the size of the reaction zone by influencing the flow field, but their direct influence on the chemical reactions is negligible. Specifically, the transport and deposition of zinc sulphide in the reactor follow the rules of a passive scalar. The mass fraction contours of zinc sulphide for the two substrates are presented in Fig. 5. The contour labels are in proportion to the magnitude of the mass fraction. For the concave substrate (left column of Fig. 5) the mass fractions close to it are small. The influence of increasing the Reynolds number is to stretch the contours in the flow direction. The consequent possibility of a higher deposition rate is, however, balanced by the prevalant stagnation region ahead of the substrate. For the convex substrate, an increase in Reynolds
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
608
0.25
0.5
1
r/D
r/D
0.5
0 0
0.5
(a)
r/D
r/D
1
x/D 0.5
1
0.25
0 0
0.5
1
0.25
0 0
1
x/D
(b)
0.5
1
x/D
(e)
0.5
0.5
0.25
10
5
r/D
r/D
0.5
(d)
0.5
(c)
3 1
0 0
1
x/D
0 0
0.25
1
0.5
0 0
1
x/D
0.25
(f )
7
1
5 3
0.5
1
x/D
Fig. 4. Contours of the production rates of zinc sulphide for the reactor with a concave substrate: (a) Re ¼ 10; (b) Re ¼ 100; (c) Re ¼ 200; convex substrate: (d) Re ¼ 10; (e) Re ¼ 100; (f) Re ¼ 200: Numbers indicate relative magnitudes.
number stretches the zinc sulphide contours towards the block (right column of Fig. 5). The absence of a stagnation region in this region suggests the possibility of higher deposition rates in this arrangement. For both geometries, Fig. 5 shows that a Reynolds number of 10 produces nearly vertical mass fraction contiurs. It is indicative of a stratification in the zinc sulphide concentration with low-concentration zones being formed in the vicinity of the substrates. From the viewpoint of uniformity of the thickness of the deposited layer, this is an advantage. The corresponding deposition rates computed as an average over the surface are expected to be low. 4.3. Wall shear stress and deposition rates It is of interest to see if the influence of the flow fields in Fig. 3 on mass transfer can be character-
ized in terms of the wall shear stresses. The shear stress distribution at the substrate surfaces are shown in Fig. 6 at a Reynolds number of 100, for two block positions. In view of the stagnation regions prevailing adjacent to the reactor tube as well as the substrate surface, these are regions of low-shear stress. At the substrate, the stress levels were seen in the computations to increase with increasing Reynolds number. At the reactor wall, the stresses were actually found to diminish with Reynolds number [18]. Fig. 6 shows that the shear stresses are higher for the convex substrate as compared to the concave. The tip of the substrate is a region of high stresses, where large velocities are encountered. A description in terms of the distribution of the wall shear stresses is useful for passive scalars such as thermal energy and non-reacting gases. Here the shear stresses can be empirically correlated by analogy to the wall fluxes of the scalar property.
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
0.25
0.5
7 9 5 10 10 3 9
5 3 5
3
1
r/D
r/D
0.5
7 9 9 5
0.25
7
0 0
1
5
3
0.5
3
5
7 1
r/D
1 r/D
0.5
1
5
3
5
9 9
7
0.25
5 1
1 0.5
(b)
0 0
1
0.5
1 x/D
0.5
3
1
0.25
9
9
5
7
1 r/D
r/D
0.5
(e)
x/D
7
9 5
0.25
3
0.5
9
0 0
1
x/D
(f )
9
5 7
3 9 1
1 0 0
1
3
3 0 0
1
x/D
9 7
3 1
7
5
0.5
7 9
0.25
1
3
(d)
x/D 0.5
0 0
1
(a)
(c)
609
3
5
7
0.5
1
x/D
Fig. 5. Contours of the mass fraction of zinc sulphide for the reactor with a concave substrate: (a) Re ¼ 10; (b) Re ¼ 100; (c) Re ¼ 200; convex substrate: (d) Re ¼ 10; (e) Re ¼ 100; (f) Re ¼ 200: Numbers indicate relative magnitudes.
The relationship between the fluxes and the stresses follows a direct proportionality rule, in the sense that regions of high shear are accompanied by high fluxes [19]. The correlation cannot be unconditionally extended to systems with reactions. The local wall fluxes of zinc sulphide are shown in Fig. 7. They are related to the material deposition rates at the curved blocks. The reactor wall and the substrate are regions of low deposition, while the edges of the substrate have high deposition rates of zinc sulphide. The wall shear stress and the wall fluxes of the transported species (and hence the deposition rates) thus show direct similarity. Hence, in the present work, the correlation between them is broadly realized under reaction conditions as well. The similarity between the wall shear stress and wall fluxes of zinc sulphide can be understood with reference to Fig. 4. Here it is revealed that the size
of the reaction zone is quite narrow and is located closer to the inflow plane of the jets. Hence, the deposition process in the reactor is dominated by scalar transport, for which the stress–flux similarity is applicable. A discussion on the relationship between shear stress and scalar wall fluxes for a reactor with a flat substrate has been presented by Saxena et al. [20]. Depositions of zinc sulphide would occur at all solid surfaces including the substrate and the tube walls. However, of interest here is the deposition at the substrate, and the effect of the surface shape on it. The local deposition rates for the convex and concave surfaces are presented in Fig. 7. It is clear that the deposition rates are higher for the convex substrate (Fig. 7b) when compared to the concave (Fig. 7a). The latter, however displays a greater uniformity in thickness, and only a marginal
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
610 0.4
0.4
0.35
0.35
concave
convex
0.3
radial coordinate
radial coordinate
0.3
0.25
0.2
0.15
0.25
0.2
0.15
0.1
0.1
Z b =0.6
0.05
Z b =0.6
0.05
Z b =1.0 o
Z b =1.0 o 0
0 0
0.05
0.1
(a)
0.15
0.2
0.25
0
0.3
0.2
0.4
0.6
(b)
Cf
0.8
1
1.2
1.4
Cf
Fig. 6. Wall shear stress profiles over the reactor substrates at Re ¼ 100; Vr ¼ 2; (a) concave, (b) convex.
0.4
0.4
0.35
0.35
0.25 0.2 0.15
Z b=1.0, Re=10 Z b=1.0, Re=100 Z b=1.0, Re=200 o Z b=0.6, Re=100
0.1 0.05 0
(a)
convex
0.3
radial coordinate
radial coordinate
concave
0
0.5
1
1.5
2
2.5
3
deposition rate (µm/min)
0.3 0.25 0.2 0.15
Z b=1.0, Re=10 Z b=1.0, Re=100 Z b=1.0, Re=200 o Z b=0.6, Re=100
0.1 0.05 0
3.5
(b)
0
0.5
1
1.5
2
2.5
3
deposition rate (µm/min)
Fig. 7. Deposition rate over the substrate as a function of Reynolds number; Vr ¼ 2; (a) concave, (b) convex.
influence of Reynolds number and the block position. For a convex substrate, the deposition rates decreases with an increase in Reynolds number. The deposition rate also decreases when the substrate is brought closer to the inflow plane (Fig. 7b, Zb =0.6). For Re ¼ 10; the deposition rate is reasonably high and the variation with the radial coordinate is uniform. The total material deposited can be calculated by averaging the deposition rate in the radial direction. The
deposited material is slightly lower at Re ¼ 10 compared to Re ¼ 100 and 200, though this is a minor disadvantage. For both geometries, the deposition rates near the edges of the substrate are high. This is because the bulk of zinc sulphide produced in the reactor bypasses the deposition surface through the annular gap. This excess material has an opportunity to stick to the corners of the substrate.
ARTICLE IN PRESS A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
In an experiemental study, Zhenyi et al. [9] have discussed the performance of a CVD reactor for zinc sulphide. The operating conditions of the reactor in the study are quite similar to the present work. In the quoted study, the deposition rate over a collection of flat surfaces oriented parallel and perpendicular to the flow is 0.3–1 mm/day, namely 0.2–0.7 mm/min. In the present numerical study, the concave substrate yielded an average deposition rate of 0.5 mm/min at a Reynolds number of 100, while it was 1.027 mm/min for the convex substrate. These values are quite close to those quoted, despite differences in the substrate geometry and orientation. The average deposition rate at a prescribed H2 S loading was computed for a convex substrate to be 1.027 mm/min compared to 0.542 mm/min for a concave substrate, an improvement of a factor of two, at a Reynolds number of 100. The doubling of the deposition rate is a significant result. The increase in the deposition rate for a convex substrate with respect to the concave is due to the following factors: (i) a better streamlining of the flow field (Fig. 3), (ii) a reduction in the size of the stagnation zone (Fig. 3), and a consequent increase in the size of the reaction zone (Figs. 4 and 5). For the concave substrate, the size of the stagnation zone barely changes with Reynolds number. Hence, the deposition rates (as seen in Fig. 7) are lower and more uniform as compared to the convex substrate, and are insensitive to Reynolds number.
5. Conclusions Fluid flow patterns and mass transfer phenomena inside a low-pressure CVD reactor have been studied via a mathematical model. Zinc vapour and hydrogen sulphide, injected along with argon react chemically in the reactor to form zinc sulphide. Of special interest is the deposition characteristics of zinc sulphide on a passive substrate held normal to the direction of the incoming flow. The following conclusions have been arrived at in the present work: 1. When a concave substrate is employed, a large stagnation region is formed where the fluid
611
velocities are small. A considerable fraction of the product species bypasses the substrate. Hence-the deposition rate is small but uniform. The deposition rate increases marginally with Reynolds number because of a slight reduction in the size of the stagnation zone. 2. A convex substrate is better streamlined, and allows the bulk of flow and the reacting species to come closer to it. This results in a smaller stagnation zone, whose size depends on Reynolds number. Species advection by the carrier gas, in addition to species diffusion significantly contributes towards the deposition of zinc sulphide on the substrate. Hence, the deposition rate increases over and above what is achieved by the concave substrate. 3. Both convex and concave substrates share unfavourable edge effects. 4. Positioning the substrate closer to the inflow plane increases the size of the stagnation zone and hence diminishes the deposition rates for both geometries.
References [1] M.L. Hitchman, K.F. Jensen, Chemical vapour Deposition, Academic Press, New york, 1993. [2] H.K. Moffat, K.F. Jensen, J. Electrochem. Soc. 135 (2) (1988) 459. [3] C.R. Klejin, Th. H. van der Meer, C.J. Hoogendoorn, J. Electrochem. Soc. 136 (11) (1989) 3423. [4] I.H. Oh, C.G. Takoudis, G.W. Neudeck, J. Electrochem. Soc. 138 (2) (1991) 554. [5] P. Duverneuil, J.P. Couderc, J. Electrochem. Soc. 139 (1) (1992) 296. [6] R.L. Mahajan, Adv. Heat Transfer 28 (1996) 339. [7] C.E. Morosanu, Thin Films by Chemical vapour Deposition, Elsevier, Amsterdam, 1990. [8] S.A. Campbell, The Science and Engineering of Microelectronic Fabrication, Oxford University Press, Oxford, UK, 1996. [9] F. Zhenyi, C. Yichao, H. Yongliang, Y. Yaoyuan, D. Yanping, Y. Zewu, T. Hongchang, X. Hongtao, W. Heming, J. Crystal Growth 237–239 (2002) 1707. [10] O. Levenspiel, Chemical Reaction Engineering, Wiley, Newyork, 1999. [11] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, Wiley, New York, 1960. [12] J.F. Thompson, Z.U.A. Warsi, C. Wayne Mastin, Numerical Grid Generation, North-Holland, Amsterdam, 1985.
ARTICLE IN PRESS 612
A. K. De et al. / Journal of Crystal Growth 267 (2004) 598–612
[13] V. Eswaran, S. Prakash, A finite volume method for Navier–Stokes equations, Proceedings of the Third Asian CFD Conference, Vol. 1, Bangalore, 1998. [14] V. Prabhakar, G. Biswas, V. Eswaran, Numer. Heat Transfer 41 (Part A) 41 (2002) 1. [15] K. Muralidhar, M. Verghese, K.M. Pillai, Numer. Heat Transfer, 23 (Part A) (1993) 99. [16] S. Kaushik, S.G. Rubin, Comput. Fluids 24 (1) (1995) 27.
[17] W.D. Seider, S.W. Churchill, AIChE J. 17 (3) (1971) 704. [18] A.K. De, Modeling of transport phenomena in a low pressure CVD reactor, Master’s Thesis, IIT Kanpur, India, 2002. [19] W.M. Kays, M.E. Crawford, Convective Heat and Mass Transfer, International Edition, McGraw-Hill, New York, 1993. [20] V. Saxena, K. Muralidhar, V. Eswaran, Sadhana 27 (6) (2002) 1.