A computational model for low-pressure diamond synthesis Citation Gavillet, Guillaume G. (1991) A computational model for low-pressure diamond synthesis. Engineer's thesis, California Institute of Technology. doi:10.7907/phqr-a780. https://resolver.caltech.edu/CaltechTHESIS:04132011-093359207 Abstract No abstract. Item Type: Thesis (Engineer's thesis) Subject Keywords: Mechanical Engineering, Engineer's thesis Degree Grantor: California Institute of Technology Division: Engineering and Applied Science Major Option: Mechanical Engineering Thesis Availability: Public (worldwide access) Research Advisor(s): Goodwin, David G. Thesis Committee: Hunt, Melany L. Leonard, Anthony Defense Date: 1 January 1991 Record Number: CaltechTHESIS:04132011-093359207 Persistent URL: https://resolver.caltech.edu/CaltechTHESIS:04132011-093359207 DOI: 10.7907/phqr-a780 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 6320 Collection: CaltechTHESIS Deposited By: Tony Diaz Deposited On: 13 Apr 2011 18:21 Last Modified: 16 Apr 2021 23:29 Thesis Files Preview PDF - Final Version See Usage Policy. 10MB Repository Staff Only: item control page

A Computational Model for Low-Pressure Diamond Synthesis Thesis presented by Guillaume G. Gavillet in partial fulfillment of the requirements for the Degree of Mechanical Engineer California Institute of Technology 199 1 11 Contents Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Index of variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Chapter 1: Introduction 7 1.1. The History of Low-Pressure Diamond Synthesis 7 1.2. Scope of this Work and Existing Research in this Area 9 Chapter 2: Methods of Low-Pressure Diamond Synthesis 11 2.1. Hot- Filament CVD . . . 11 2.2. Combustion Torch CVD 12 2.3. Plasma CVD 13 Chapter 3: Analytical Formulation 15 3.1. Similarity Solution for the Flow Field 15 3.2. Boundary Conditions. . . . . . . . . · 20 3.2.1. Boundary Conditions for the Flow Field · 20 3.2.2. Inlet Conditions for Chemical Species · 21 3.2.3. Wall Condit ions for Chemical Species · 21 3.2.4. Boundary Condition for Wall Axial Velocity 3.3. The Radial Pressure Gradient Parameter 22 23 Chapter 4: The REAC5 code for CVD simulation · 25 4.1. Method of solution .. . • _C) ?- 4.1.1. General Principle · 25 4.1.2. Two-Part Solution Technique · 26 III 4.2. Damped Newton method. · 27 4.3. Time relaxation · 30 4.3.1. Transformation to an Unsteady Problem · 30 4.3 .2. The Problem of Stiffness · 30 4.3.3. Modified Boundary Condition for the Mass Fractions. 32 4.3.4. Modified Boundary Condition for Wall Axial Velocity 33 4.4. Ini tial Guess for the Newton Method . 33 4.4.1. Flow Field . . . . . . . . . . 34 4.4.2. Pressure Gradient Parameter 34 4.4.3. Chemical Species . 35 4.5. The Computational Mesh 36 4.5.1. Fixed nonuniform grid . 36 4.5.2. Adaptive gridding . . . . 38 4.6. Problems Specific to the Last Node . 40 Chapter 5: Results .. 5.1. Validation of the Code · 42 · 42 5.1.1. Flow Field Calculation · 42 5.1.2. Species Profiles Calculation · 43 5.2. Hot Filament CVD Simulation. . .44 5.2. 1. Description of First Experiment .45 5.2.2. Description of Second Experiment . · 48 5.2.3. Chemical Mechanism for Simulation · 48 5.2.4. Results. . . . . . . . . . . . . .. . · 49 5.3. Preliminary Plasma Torch CVD Simulation · 53 5.3.1. Experiment · 53 5.3.2. Chemical Mechanism for Simulation · 54 IV 5.3.3. Results . . . . . . . . . .. . . . . . . . . . . .. .. . . . . 54 Chapter 6: Furt her Improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Chapter 7: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Appendix A: Mass Flux Balance at the Substrate Surface .. . . . . . . . . . . . . . . . . 63 Appendix B: Chemical Mech anism File . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 -1- 2 Acknow le dgements I would first like to thank my adviser, Dr. David Goodwin , for his useful and extensive advice, which helped me solve the multiple problems associated with t he completion of this work. Thanks also to Dr. Melany Hunt, who allowed me access to the THERMI node, sparing me many long waits for resul ts on slower systems or costly CADRE bills. Prof. Anthony Leonard was of considerable help in the solution of numerical problems which arose during the development of the code. J oe Sivo gave me useful hints for t he analytical formul at ion of the problem. Angela Shih h elped me master the AUTOCAD software, with which some of t he figures in this thesis have been drawn . I finally wish to thank my office-mates, Susan Beatty. Ann Chopra, Dave Kaufman , Sanjay Kumar and Susan Melnik, who contri buted to making my workplace a very pleasant one. Guillaume Gavillet 3 Index of variables ROlnan: a parameter for axial coordinate transformation bj j th coefficient of interpolation parabola Ck mean molecular velocity of species k Cpk specific heat at constant pressure of species k Cp mean specific heat at constant pressure Dk mass diffusi vi ty of species k E partial equilibrium ratio for hydrogen recombinat ion (section 5.3 .3) hk specific enthalpy of species k H height of CVD reactor kB Boltzmann constant, 1.381 X 10- 16 erg/ K mk mass of an atom of species k lvh molecular mass of species k M mean molecular mass AI! a Mach number nk number density of species k NA Avogadro's number, 6.02252 X 1023 Ne/ number of elements in the chemical mechanism Neq number of equations to be solved per node Nm number of mesh points in the computational domain Nspec number of species in the chemical mechanism Nu number of unknowns of the problem p pressure P component of pressure depending on z only Pe P eciet number Pr Prandtl number r radial coordinate R gas constant (average density) Re R eynolds number 4 S scaling factor for chemistry time scale (section 4.3.2 ) sl standard liters (1 atm and 298 K) t time T temperature u axial veloci ty Uk diffusion velocity of species k v radial velocity V l/r times radial velocity w circumferential velocity W angular velocity x transformed axial coordinate X k mole fraction of species k Yk mass fraction of species k Z axial coordinate Zj axial coordinate of node i Calligr a phic : A absolute component of error tolerance Ck surface creation rate of species k 'Dk surface destruction rate of species k :T component of Jacobian matrix Nk number of hydrogen atoms in the species k n relati ve component of error tolerance 5 Greek: 'Yk surface reaction probabili ty of species k r ratio of specific heats o tolerance for grid ada ptat ion c: tolerance for convergence C. Zi+! - Zi K: mean thermal conductivity ,\ damping parameter for Newton method A l / r times r adial pressure gradient, assumed constant J1.k dynamic viscosity of species k J1. mean dynamic viscosity Wk molar product ion rate of species k n ¢ substrate angular velocity ¢i ith unknown Wk molecular flux of species k p density () circumferential coordinate vector of unknowns Elk t hermal diffusion ratio of species k Xj per t urbation for numerical differentiation of the variable j Y van der Monde determinant 6 Subscripts: cal c calculated ch characteristic of chemistry variables eq equilibrium JI characteristic of flow-field variables int intermediate L at inlet res residence ( time) speciJ specified W at wall W + 1 next to wall o steady-state 6. relati ve to convergence of Newton method X relative to numerical differentiation Superscripts: + with positive z-component of molecular speed with negative z-component of molecular speed ( n ) at iteration n s surface 7 Chapter 1. Introduction Diamond has long been known to be an almost ideal material for a wide- ranging variety of applications, from wear-resistant cutting tools to semiconductor circui ts. It has a high thermal conductivity (up to 20 W /cm/K, at least twice as much as copper or silver), it is the hardest known material, it is optically transparent over a broad range of wavelengths, it has a fairly high breakdown voltage (70000 V j cm ) and can be doped to form a semiconductor (although only p-type doping has been achieved to date). The only factor currently limiting its use is its high price due to its rarity. The known supplies of natural diamond do not suffice to satisfy the ever-growing needs of the industry, and the current high-pressure techniques for artificial diamond synthesis, although they yield gems of remarkable purity, can barely compete wi th extraction of natural diamond as far as price is concerned. The range of possible industrial applications is impressive, and could lead to major breakthroughs. For example, a diamond-coated tool would offer a much higher resistance to wear, and therefore a h igher lifetime, than the hardest current tungsten carbide tools . A diamond coating on an optical window would enhance its resistance to scratching without appreciably altering its transparency. A bulk diamond lens would have a much better heat dissipation rate than an ordinary lens, therefore permitting easier use in powerful lasers. Cheaper heat sinks could be manufactured for the semiconductor industry. If the technique of hetero epita:{ial growth is mastered, the way would be open to the production of diamond-based semiconductors. All these applications (except heat sinks) require the production of a thin film rather than bulk diamond. Therefore, it is not surprising to see some of the major industrial powers-Japan, the United States, the Soviet Union- fund significant research programs in diamond thin films. Introduction 8 1.1. The History of Low-Pressure Diamond Synthesis The possibility of diamond synthesis was suggested for the first time in the beginning of the twentieth century. It was initially believed that high pressures and temperatures were necessary, since the diamond form of carbon was thermodynamically stable only in this region. Increasing knowledge of the phase diagram of carbon, due to the work of Leipunskiyl among others, helped define the range of diamond stability, but also showed the existence of a domain of metastable diamond existence near or below atmospheric pressure: This led to the thought that • low-pressure diamond synthesis was also possible. Research in artificial diamond growth began in Sweden! and in the United States shortly after the end of World War II. It is widely believed that the first efforts to achieve controlled diamond growth and production were all focused on high-pressure techniques. However, it was the research group of Eversole,6 which was the first to succeed in the beginning of 1953 with a low-pressure process, a few months before a Swedish team at Allemanna Svenska Elektriska! reported the first high-pressure diamond synthesis. This illustrates the fact that research in low-pressure diamond growth is far from being a new topic, even if, unlike its high-pressure cousin, it is just entering the industrial stage. Eversole succeeded in depositing a tilln film of diamond and graphite on a diamond substrate by exposing it to a Illgh-temperature, low-pressure flow of carbon monoxide or "any gas decomposable to a methyl radical,,6 -with hindsight, a rather prophetic statement. To remove the undesired graphite, he had to periodically expose the film to a high-temperature ( 1300 K), high-pressure (50 bar) hydrogen environment. Tills cyclic process yielded films of acceptable purity, but the growth rate was far too low (0.01 J.lm/hr) for any industrial application to be considered. Moreover, he did not report growth on a non-diamond substrate, which rendered the interest in his method purely theoretical. Research in this area was discontinued in the United States shortly afterwards. Introduction 9 A Soviet research team under Derjaguin,1 working independently from Eversole's group, also reported low-pressure diamond growth on a diamond substrate as soon as 1956. For eighteen years, the Soviet Union remained the only country with a known interest in low-pressure diamond synthesis, but little progress was reported. The real effort in this field started in 1974 with the launch of a large-scale research program by Japan's National Institute for Research in Inorganic Materials (N IRlM) . Derjaguin's team reported the first growth on a non-diamond substrate in 1981, but it is the Japanese who made the most significant contributions, bringing low-pressure diamond synthesis to the industrial stage. Matsumoto et al. 8 developed the hot-filament chemical vapor deposition (CV D), then the plasma CVD processes. Growth rates of plasma CVD at atmospheric pressure have recently surpassed 500 J.lm/hr.9 Research in Japan also resulted in the discovery of the combustion torch deposit ion process. The first industrial applications-diamond-coated tweeters for sound systems, diamond loudspeaker cups, X-ray windows-were presented in 1988. The market for diamond thin film technology is expected to grow very quickly in the next decade, with an estimated! $885 million worth of sales by the year 2000. 1.2. Scope of this Work a nd Existing R e search in this A r ea Understanding diamond CVD requires a good knowledge of the flow field m the reactor and of the mechanism of transport of carbon atoms tp the substrate. Velocity measurements in a CVD reactor are difficult, and the chemical mechanisms involved are not yet fully understood. This prompts the development of computational models to test the validity of new hypotheses on the flow field and on the chemistry against experimental results . Unfortunately, in the present state of computer technology, the full simulation of CVD in an actual reactor, where the flow is three-dimensional (sometimes recirculating) and whose chemical mechanism involves a large number of species, is of prohibitive cost. Therefore, simplifications have to be made on reactor geometry, chemical mechanism and flow field properties. Introduction 10 The known CVD codes are not numerous; all include such simplifications, and can be divided into two groups. A first approach is to concentrate on the flow field properties. In a typical CVD experiment , a small amount of reactants is diluted into a carrier gas, or one reactant is present in overwhelmingly large quantity, which means that the transport properties of the mixture can be approximated by those of the dominant specIes. Moreover, the heat release due to the reactions has little influence on the temperature distribution in the reactor. Therefore, neglecting the chemical mechanism altogether and solving for the flow field with the pure carrier gas gives a good estimate of the actual conditions. Starting from this hypothesis, Evans and Greifl2 developed a two-dimensional axisymmetric code, including finite-radius effects and buoyancy forces, which accurately calculates velocities and temperatures in a finite-radius reactor. Another possibility is to concentrate on an accurate description of the chemical mechanism while simplifying the fluid mechanics, transforming the 3D problem into a 2D or ID axisymmetric problem by neglecting, for example, buoyancy and finiteradius effects. If needed, a physical experiment can be devised whose geometry is close to these ideal conditions. This is the approach of Coltrin et al.,ll who developed a one-dimensional axisymmetric laminar code including gas-phase and surface reactions, but only applied it to silicon CVD. A comparison with results from E vans and Greif in the case of an aq;on flow over a rotating substrate showed that the ID approximation is valid over 90% of the subst rate radius. The scope of this report is the development of a cost-efficient code applicable to the simulation of diamond CVD in any axisymmetric reactor whose height is small compared to its diameter. This means that an approach similar to Coltrin's has been chosen, with emphasis being put on chemist ry rather than on the fluid flow. The code is applicable to many different reactors, for example hot-filament , flatflame and plasma reactors (to be presented in more detail below), and therefore, we will focus on these methods. The principle of the code will be discussed , then results of flow-field validation calculations and comparisons of diamond CV D simulations with experimental data will be presented. Chapter 2. Methods of Low-Pressure Diamond Synthesis The known methods of low-pressure diamond synthesis (see Spear,2 Angus & Hayman 3 for a review) are numerous. We will focus on three main techniques: hotfilament, combustion torch and RF plasma torch, which will be briefly described. 2.1. Hot-Filament CVD 1_ Hz , CH. 1 300 K --+ 200 \,/ / / ....-- / filo.Ment \ 2500 K "- "- ' - ~ substrate 1200 K Figure 2.1.1: Hot-filament CVD. The principle of the hot-filament CVD method (HFCVD ) is to thermally decompose a carbon-containing gas to deposit solid carbon on a substrate. More specifically, a mixture of hydrogen and methane, in the typical (molar) proportions of more than 99% H2 , is injected in a low-pressure chamber ( typically 20 to 100 Torr) and flows past a heated filament (usually tungsten, at a temperature of 2200 Methods of Low-Pressure Diamond Synthesis 12 K or more) before impinging on a substrate (silicon, quartz, molybdenum , platinum, diamond ... ) a short distance (typically 5 to 15 mm) away from the filament (see Fig. 2.1.1 ). The exposure to the filament triggers a pyrolysis mechanism as well as partial dissociation of molecular hydrogen , which leads to the formation of various hydrocarbons. Contact with the substrate surface results in the decomposition of some of these species, whose carbon atoms are deposited on the surface either as graphite, amorphous carbon or diamond. However, it has been shown that graphite is etched away faster than diamond when in presence of significant amounts of atomic hydrogen (this fact was used by Eversole in the first diamond CVD reactor). In a hot-filament reactor, dissociation of molecular hydrogen occurs on a fairly large scale near the filament, and the substrate is still hot enough (typically 1000 to 1200 K ) to prevent H2 recombination. Therefore, it is not surprising that diamond films grown in a HFCVD reactor can be of high quality, with Raman spect ra almost as good as natural diamond's . This is counterbalanced by the rather low growth rates (typically 1 to 10 Jim/hr) allowed by this method. Adding small amount s of an oxygen-containing compound (carbon monoxide, water, atomic oxygen) can boost growth rates to 40 Jim/hr 1 with only a small penalty in purity. 2 . 2. Comb u stion Tor ch CVD Combustion torch chemical vapor deposition (CTCVD) does not reqwre an external heat source like HFCVD . The react ion mechanism involved is no longer a pyrolysis, but a true combustion mechanism. A hydrocarbon (often acetylene) is burned in oxygen, the flame being run slightly fuel-rich (fuel-air equivalence ratio between 1 and 1.2). A nondiamond substrate is held in the "acetylene feather, " the intermediate region of the flame where the excess acetylene reacts with air which diffuses from outside the flame (see Fig. 2.2.1 ). The substrate has to be water-cooled in order to keep its temperature in the optimal range. This method works at low or atmospheric pressures (where it is most efficient). Therefore, nothing more than an ordinary welding torch is required to generate a diamond film at a higher growth rate than with HFCVD (up to 60 Jim/hr4 ), but on a smaller area. The potential Methods of Low-Pressure Diamond Synthesis acetylene feather 13 s toichioMe tric core substrate 1200 K Figure 2 _2_1: Combustion torch CVD. of this method is considerable: progress in growth rate, film purity and operating temperatures would make it possible for a worker to perform diamond coating of a given surface simply by exposing it to his welding kit for a certain time period. A variation on this technique is now under investigation, in which a low-pressure flat-flame burner is used instead of a welding torch. 2.3. Plasma CVD It seems natural to think that a plasma, in which dissociation occurs on a large scale, should be an environment favorable to diamond deposition. A gas flow (hydrogen, a hydrocarbon and sometimes a carrier gas), heated to plasma temperatures (typically 4500 K) by exposure to radio frequency waves, microwaves or a DC current, impinges on a water-cooled substrate, creating a thermal boundary layer around it (see Fig. 2.3.1 ). Cooling occurs fast enough for atomic species to exist in super-equilibrium concentrations, and favorable conditions for diamond growth Methods of Low-Pressure Diamond Sy nthesis 14 o o o RF coil o o o aio.MOna ~ substro.te 1400 K Figure 2_3_1: Plasma torch CVD. exist. P lasma methods work well at atmospheric pressure, and yield by far the most impressive growth rates (more than 500 p.m/hr9 with a DC plasma torch). 15 Chapter 3. Analytical Formulation 3.1. Similarity Solution for the Flow Field The theoretical formulation of any of the three CVD techniques described in Chapter 2 is quite similar. First of all, we are dealing with a steady flow . In all t hree cases (HFCVD , CTCVD, RFPCVD ), especially in the case of flat-flame CTCVD , the substrate diameter and the gas inj ector dimensions are large compared to t he distance H between filament (HFCVD ) or gas inlet (CTCVD , RFPCVD ) and substrate. Therefore, we will consider the substrate radius to be infinite with respect to H . Furthermore, we will assume that t he gas is injected at a constant velocity UL, uniform throughout the injection surface. 'vVe also consider an axisymmetric reactor, i.e. , a disc substrate below a circular injection orifice . The geometry of the idealized reactor, along with t he cylindrical coordinate system used and the definit ion of the velocity components, is shown on Fig. 3.1.1. I [J_[IJ_[IJ_[IJ_[j~ t I I I HI I : ~ _ ru 01 ®w.e ~ ,r T" Figure 3.1.1: Idealized reactor, coordinate system and notations. 16 A nalytical Formulation The substrate can be fixed or rotating with a uniform angular velocity Sl. It has been observed in silicon CVD experiments 11 that the thickness of the deposi ted film is more uniform when the substrate is rotating. Moreover, we will neglect buoyancy effects. In all CVD conditions for which simulations were performed , the height-based Reynolds number: was lower than 400. Therefore, the flow is laminar. We solve for the mass fractions Yk of all N. pec species but one (since t he sum of all mass fractions is equal to 1, there are only N. pec -1 independent mass fr actions), for the velocity components v, 1.1, w, and for the density p, since it is a compressible flow. In a system of cylindrical coordinates (r, B, z) as shown on Fig. 3.1.1 , and with the axisymmetric hypothesis, the conservation equations are: Conservation of species k (k = 1, N. pec - 1): (3.1.1 ) Radial momentum: ov w2) =--+op a [ ov 2 (ov v -au)] p(v-ov -or +-+ or +u--oz r or or 2Jl-or +-Jl 3 r oz +~oz [Jl (auor +ov)] oz +2Jlr (ov or _~) r (3.1.2) Circumferential momentum: p( v-OWor +u-owoz +-Vw) ~or [Jlr~or (W)] + ~oz (Jl OWoz ) + 2Jl~or (W) r = r r (3.1.3) Axial momentum: p op +ora [(au (vau or + au) oz =- oz Jl or +Ov)] oz +2Jl-oza [au -oz - -31 (ov -or + -vr + -au)] oz +-Jlr (au -or + -ov) oz 1.1 (3.1.4) A nalytica l Formulation 17 Enthalpy: p L ah) = v-ap + -a (aT) P ( v-ah + uK;+ J.I (av)2 - - a ( p N. " Uk h kY k) ar az ar az az az az (3.1.5) k=l Continuity: a(pv ) + a(pu) + pv = 0 ar az r (3.1.6) The velocity Uk appearing in Eq. (3.1.1 ) is the species cliffusion velocity, incorporating both mass and thermal cliffusion (Soret effect): ~~ UkYkM [Mk0 k a (lnT ) _ aCMYk ) ] a a z z (3.1.7) This is only an approximation, since mixture-averaged mass diffusivities are used instead of binary cliffusion coefficients. We also need the equation of state: p=pRT (3.1.8) where R is the gas constant. With the assumptions mentioned at the beginning of this sect ion , it is possible to postulate a similarity solution for the raclial and circumferential components of the velocity: w(r,z) = rW (z) v(r,z)=rV(z), (3.1.9) We will therefore solve for V and W instead of v and w. Let us make the adcli tional hypotheses au = 0 ar ' aT =0 ar ' - aYk = 0 ar and consider the enthalpy equation (3.1.5). Let us recall that N,pec h= L hkYk k=l (3.1.10) Analytical Formulation 18 and that , therefore, oh oz Also Eq. (3.1.5) can then be rewritten: Taking into account the species conservation equation (3.1.1 ), this becomes: The right-hand side terms can be neglected with respect to the convective term. Let us compare the pressure gradient term with the convective term. Using the approximations op - oz u}, ~p ­ 2H !!.- ("KdT) ~ "KTL dz dz H2 we obtain the following condition for the pressure gradient term to be small: Analytical Formulati on 19 which can be rewritten as: and, since this yields: 2 1 M a «: (r _ 1) Pr which holds in all the cases studied. Therefore, we can neglect this term. Similarly, to neglect the velocity gradient term , we need: 2 M a 2 «: - - - - - (r -l ) RePr which holds in all the cases studied, and we can neglect the term. Finally, using the hypotheses (3.1.9) and (3.1.10), and neglecting the right-hand side of Eq. (3.1.11 ), the system (3.1.1 ) to (3.1.6) becomes: Conservation of species k: (3.1.12) Radial momentum: pu -dV + p (V 2 - W 2 ) - -d ( f.J. -dV) + -1 -op = 0 dz dz dz r or (3.1.13) Axial momentum: (3.1.14) Circumferential momentum: d ( f.J. -dW) =0 pu -dW +2pVW-dz dz dz (3.1.15) Thermal energy: (3.1.16) Continuity: du udp - +2V+ - - =0 dz p dz (3.1.17) Analytical Fonnulation 20 This system can be rendered fully one-dimensional , depending only on the axial variable z, if we make the additional assumption lap -a = A = constant r r or, equivalently: r2 p(r,z) =P (z)+ A 2 (3.1.18) If, moreover, P( z) :::p Ar 2 /2 (confirmed a posteriori), then the equation of state can be rewritten as: P (z) = pRT(z) (3.1.19) With a few simplifying hypotheses, we have therefore transformed the axisymmetric two-dimensional problem of the flow inside a CVD reactor into a onedimensional problem described by the system (3.1.12) to (3.1.19), which can be solved more easily. 3 .2 . Boundary Conditions 3.2 .1 Boundary Conditions for the Flow Field To complete the formulation of the problem, it is necessary to state boundary condi tions for flow-field and species equations. The flow-field conditions for a reacting flow are not absolutely identical to the ones for a nonreacting flow. Indeed, since we are considering deposition on the substrate surface, there must be a nonzero normal flux of some species, and t herefore a non-zero substrate axial velocity uw. The condition on uw will be t reated in section 3.2.4. The other boundary conditions for flow-field variables are the usual ones: Analytical Fonnulation 21 at the substrate: V (z = 0) = 0 W(z=O)=Sl T(z=O)=Tw at the inlet : V (z = L ) = 0 u(Z=L)=UL W(z=L)=O T (z =L )= TL P(Z=L)=PL All these quantities are problem-specific, and are set by the user at the beginning of the calculation. The condition on the inlet pressure serves as boundary condi tion for the continuity equation (3 .1.17) after t ransformation by the equation of state (3 .1.19). 3.2.2 Inlet Conditions for Chemical Species The composition of the gas feedstream at the reactor inlet is fixed. However, it is not necessarily the equilibrium composition at the inlet temperature. For example, in the case of HFCVD, the gas passing by the filament does not have time to reach thermal equilibrium, and enters the reactor at a lower temperature which will be chosen as the inlet temperature TL. Moreover, the residence time of the gas at this temperature is too short for the gas to reach chemical equilibrium. The inlet temperature TL is determined by experimental observation, and the residence time is chosen so as to yield species concentrations close to measured values. A preliminary calculation is performed by integrating the rate equations: [} [Xkl at = Wk (3.2.1) over the residence time, keeping the temperature fixed at TL. The inlet mass fractions calculated this way are then kept constant throughout the whole computation: Yk ,L = Yk (TL , tr •• ,d 22 Analytical Formulation 3.2.3 Wall Conditions for Chemical Species An accurate description of the diamond CVD process would require detailed knowledge of the surface chemistry on the substrate, which is not available at this time. It is not even known with certainty which species are precursors to diamond deposition. A first crude attempt to ascertain the role of various species is to introduce the concept of "surface reaction probability." We will assume that a molecule of a carbon-bearing species k impinging on the substrate surface has a probabilitY'Yk of reacting, integrating all its carbon atoms into the crystalline lattice and releasing its hydrogen in molecular form. A mass flux balance (see Appendix A) yields the following condition: Mk kBTw C. ( 'Yk) ] ['Yk-Ck,W 4- + 1 -"2 Uw Yk,W M Pw _ (1- ~) D;: k [Mkek(d(~:T))w - e~Yk) )J (3.2 .2) where and q is the surface molar creation rate of the species k (zero for all species except H2)' This is the wall boundary condition for the species k . 3.2.4 Boundary Condition for Wall Axial Velocity Equation (3 .2.2) can be rewritten in the form: (3.2.3) Summing over all species and making use of the property N,pee L YkUk = 0 k=l Analytical Fonnulalion 23 yields an expression for the wall bulk velocity resulting from the deposition process: PWUW = N. p " C'M v k£.. ' " k k - PW' k ,Wi 4 L.. 1 _ li k=l 2 (3.2.4) which will be used as boundary condition for the wall axial velocity. We now have a system of Neq = N,pee + 4 equations (N. pee - 1 species conservation equations and 5 flow-field equations), and the correct number (2Neq - 1) of boundary conditions with respect to the order of the equations. Yet the problem is not completely formulated, since nothing has been said of the pressure gradient parameter A appearing in Eq. (3.1.18). 3.3. The Radial Pressure Gradient Parameter In case of a nonrotating substrate, the existence of a favorable (i.e. , negative in the coordinate system of Fig. 3.1.1 ) radial pressure gradient is a necessary condition for a flow to occur. In the case of a rotating substrate, however, it is also possible for a flow to occur in the absence of a radial pressure gradient (due to viscous entrainment by the substrate, as in the von Karrruin l 6 solution for a viscous stagnation flow on an infinite rotating plate in a semi-infinite domain) or even with an adverse (positive) gradient. With the hypothesis (3. 1.18) on the static pressure in the reactor , it is the parameter A that determines the sign of the radial pressure gradient. It is clear that for a given substrate rotation rate n (including the nonrotating case n = 0), there is a one-to-one relationship between the radial pressure gradient and the quantity of fluid "pumped" by the plate (n # 0) or drawn into the reactor (n = 0). In other words , the value of A and of the inlet axial velocity UL are interdependent. When viscous entrainment is present, it is possible to define an "ideal" inlet axial velocity corresponding to a zero radial pressure gradient, and it is clear that any inlet velocity superior (in modulus) to this ideal value will yield a "forced" flow (negative gradient), whereas any inferior value (in modulus) will yield a "starved" flow (positive gradient). Analytical Fonnulaiion 24 Therefore, the parameter A is an eigenvalue of the problem: for a given set of boundary conditions, only one value of A exists that yields a solution. Conversely, if the boundary condition on the inlet axial velocity tiL were replaced by a condition on A, tiL would then be uniquely determined. 2S Chapter 4. The REAC5 code for CVD simulation 4.1. Method of solution 4.1.1 G e neral Principle We now want to solve the previously described system of Neq orclinary clifferential equations in the computational domain [z = 0, z = H ]. The first step is the cliscretization of the domain into a finite number N m of points or "nodes." Then, the equations are finite-clifferenced at these points. The new system to be solved is a set of N. q x N m = Nu algebraic equations. As we saw, the parameter A appearing in Eq. (3 .1.13) is an eigenvalue of the problem, and will be calculated iteratively. For this purpose, the boundary condition u(z = H ) = UL ,speci! is abandoned. An initial value of A is picked, which yields a flow field with an inlet axial velocity UL ,calc ' Then A is iteratively adjusted (multiplied by the ratio UL,speciduL ,calc ) until (usually E = 10- 2 ) . It can happen that the actual value of A is bracketed by two values AI , A2 at the last two iterations (with corresponcling inlet axial velocities UL I, uL2, the index 2 indicating the last iteration performed), and that the correction performed on Al yields A2 and vice verJa. If this happens , then a special correction technique is applied: A2 is multiplied by the ratio which yields an intermediate value of A and makes convergence possible. The system will be solved by a combined damped Newton iteration-time relaxation method similar to the technique used in the one-dimensional flame model developed by Grear et al. lo and in the axisymmetric CVD reactor model developed The REAC5 code for CVD simulation 26 by Coltrin et a/. l l An initial guess is made on the unknowns , and an attempt is made to solve the (steady-state) system by the damped Newton method . If this does not succeed, the problem is transfonned into an unsteady one by adding the relevant time derivative to the equations of conservat ion (3.1.12) to (3.1.17), and by calling a time integration routine for a fixed duration. When time integration is completed, the system is converted back into a steady one and a new attempt is made with the Newton method. This combination has the advantages of both techniques ( the high speed of the Newton method and the loose restrictions on the initial values of time integrat ion), and minimizes the overall convergence time. 4 .1.2 Two-Part Solution Technique A distinction should be made between "flow field" equations (3.1.13), (3.1.14), (3.1.15), (3.1.16), (3.1.17) , and species conservation equations of the type (3.1.12) . Since the composition of the gas in the reactor is usually dominated by an overwhelmingly present species (either molecular hydrogen or a "carrier gas" like argon), the quantities of reactants are comparatively small. Therefore, the heat release or absorption due to chemical reactions and the change in the gas mean transport properties will not have a major influence, and we can expect density, velocity and temperature profiles in the reactor to be very close to those obtained with the pure, non-reacting carrier gas or with "frozen" chemistry (i.e., performing a nonreacting calculation with a constant chemical composition throughout the computational interval). This argument has been used by Evans and Greif,12 who only considered the carrier gas in their fully three-dimensional calculation of pressure, velocity and temperature inside a finite-radius CVD reactor. This property is of importance for our purposes. Indeed, solving for the flow field only would require little CPU expense, and since the flow field is nearly unchanged when reactions are included, the result of this nonreacting calculation would be an excellent initial guess, obtained at low cost, for the full (reacting) calculation. This should make convergence of the Newton method significantly faster . Therefore, we split the computation in two phases. First, chemical reactions are turned off by replacing the species conservation equations (3.1.12) by Yk = constant, 27 The REAC5 code for CVD simulation and a solution is obtained for velocity, density and temperature fields and for the parameter A. Then reactions are allowed to proceed by using Eq. (3.1.12 ) again , and the final solution is obtained. During the development of the code, it has been observed that convergence of the reacting calculation is slower when thermal diffusion is included in Eq. (3.1.7). Also, the reacting calculation almost always includes two iterations on the pressure gradient parameter A. Turning off thermal diffusion in the first reacting iteration on A, by set ting all thermal diffusion ratios to zero, has proven to be a good way of achieving quicker convergence while taking thermal diffusion into account. The thermal diffusion ratios are restored to their normal values at the beginning of the second reacting iteration on A. If convergence of the reacting calculation is reached after only one iteration, a second iteration including thermal diffusion is forced. 4.2. Damped Newton method The Nu unknowns are arrayed into a single vector: ;j = (YI ,I, ... , YI ,N",, -b VI, UI, WI, T I , PI, Y2 ,1,· ··,PNm -l, YNm,l, ... , YNm ,N"pee-1, VN m ) UNm ) WN m ) TNm ,PNm) The system is now of the form F(;j) = 0, adequate for the use of a Newton algorithm. Let ;j(n) represent the approximate solution of the system after n iterations. In the "classical" (i.e., undamped) Newton method, the next approximate solution ;j(n+l) is obtained by: (4.2.1 ) where [J'(n)j is the Jacobian matrix: An initial guess ;j(0) is necessary to start the computation. If it is too far from the solution, the residuals F( ;j(0» will be very large, and it is likely that the 28 The REAC5 code f or CVD simulation Newton iteration will yield a vector i even further away from the actual solution. There exists a "domain of convergence" of the Newton method, analogous to the domain of stability of integration schemes, in which the initial guess should be made. The mathematical determination of this domain is a formidable task except for very simple problems, and therefore the chosen initial values of the unknowns are usually the result of an educated guess. In our case, it is possiblel 5 to devise a simplified flow model which yields a good initial guess for all CVD experiments. This will be discussed in section 4.4. It can happen that an iteration of the form (4.2.1 ) hinders convergence, i.e., the norm of the F-vector increases after the (n + 1)th iteration. In this case, it is often helpful to "damp" the iteration by reducing the increment on each unknown 4>i. This can easily be done by introducing the "damping parameter" A (0 < A ::; 1), and modifying Eq. (4. 2.1): ;j<n+l) = ;j<n) _ A [.J(n)] - I F (;j<n)) = ;j<n) _ At:.;j<n) (4.2.2) The evaluation of the Jacobian matrix [.J(n)j is the most time-consuming part of the computation. Fortunately, it is not necessary to go through the computation of all coefficients. First, let us consider one line of the Jacobian. It contains the partial derivatives of one equation at one node with respect to all unknowns at all nodes. If this line is, say, for the node i, only the unknowns at the nodes i-I, i, i + 1 may appear in this equation due to the nature of the finite-difference schemes used (see section 4.5). Therefore, the Jacobian matrix has a banded structure. Of all 3Neq variables in the band, only a few appear indirectly (e.g. temperature-dependent transport properties) or directly in the equation (see for example Eq. (3. 1.13)), which means that most of the coefficients inside the band will be zero . A large amount of CPU time can be saved by recognizing in advance which coefficients will be non-zero and by computing only these. This approach is used in the code, and the non-zero coefficients are computed numerically by the usual differentiation formula: 'T. _ vi, - Fi (t/>J, ... , 4>j- l , 4>j + Xj, 4>j+l, ... , 4>NJ - Fi (4) 1, ... , 4>j-l , 4>j, 4>j+l, ... , 4>NJ Xj 29 The REAC5 code for CVD simulation where Xj = AX + Rx i¢>ji, Ax and RX being constants set by the user at the beginning of the calculation. All computations presented here were performed with Ax = 1.01 X 10- 7 and RX = 9.9 X 10- 9 . The value of Ax could seem too high in regard of the low mass fractions of some species, but it should be remembered that a well-known feature of Newton methods is that they work well even if the Jacobian is not accurately computed (however, zeroes and non-zeroes have to be accurately positioned). Due to this resistance of Newton methods to Jacobian inaccuracies, it is not necessary to compute the Jacobian at each iteration. Instead, it is computed once at the beginning (in our notation, this is [.7(0)]), then updated only when needed , i.e., when convergence can no longer be obtained with the old matrix. The acceptance of the solution vector at the nth iteration depends on a "look ahead" procedure. The new solution vector ;j(n+l) is first computed via Eq. (4.2.2) with>' = 1, then a tentative 6.;j(n+l ) is computed. The iteration is deemed complete if max 15;i5;N" 6.¢>(n+l) IA,; + R,;i rjJ~n+ l )i S I< max 15;i5;N" I 6.rjJ(n) t A,; + R,; i¢>~n)i I ( 4.2.3 ) where A,; and R,; are constants chosen by the user at the beginning of the calculation. All computations presented here were performed with A,; = 10- 19 and R,; = 10- 4 If the condition (4.2.3) is not satisfied, the damping parameter>' is halved. This allows j(n+l ) to be closer to j(n) , which in turn reduces 6.j(n+I). The reduction of >. continues until either (4.2.3) is satisfied or >. falls below a specified value (in our case, 2- 10 ) . If convergence has not been obtained after 10 reductions of >., a new Jacobian matrix is computed and the iteration is restarted with >. = 1. If convergence still cannot be obtained, time relaxation is called. It has been observed 10 that the first Newton iteration after a Jacobian update nearly always yields a j which is "further away" from the solution (bigger residuals) than the previous one. Yet during subsequent iterations, the residuals may decrease again. Therefore, the "look ahead" criterion has to be bypassed during the first iteration in order to avoid an unwanted rejection. 30 The REAC5 code for CVD simulation 4.3. Time relaxation 4.3.1 Transformation to an Unsteady Problem Let us consider a well·posed steady-state problem i"{ ¢) - 0 with a umque solution ¢o. Consider now the associated time-dependent problem: -8¢ =F~ (~) ¢> &t If the problem is such that ¢ approaches an equilibrium solution ¢eq as t ----+ 00 and that the order of the system is not modified when adding the time derivative, then ¢eq has to be equal to the steady-state solution ¢o, since the time-dependent problem to be solved approaches the steady-state problem more and more closely (uniform convergence). Hence, it is possible to solve the steady-state problem by transforming it into a time-dependent one and integrating until equilibrium is reached. 4.3.2 The Problem of Stiffness This approach looks much more tolerant than a Tewton method since there is no domain of convergence in which to start to obtain a solution. However, stability requirements are a limiting factor in the choice of the time step. Within these boundaries, the time step can usually be adjusted to maintain the local error on the solution within a tolerance specified by the user. Unfortunately, when the unknowns vary over dramatically different time scales (a typical feature of most chemical kinetics problems), some of the eigenvalues of the Jacobian matrix can have very large moduli, which makes the stability requirements of classical integration schemes, like Adams-Bashforth or Runge-Kutta, too stringent for time step adjustment routines to be efficient. This problem of large eigenvalues is called stiffness, and it affects nearly all chemical kinetics simulations. Moreover, since a high time derivative is a source of error, the time step should be kept minimal when any of the components of the system varies quickly. In a problem where sharp variations occur at different times, this means that the time step, already severely limited for stability reasons, will have to be kept at a very low value during most of the computation, The REAC5 cod, for CVD simulation 31 thus making integration long and costly. There exist methods specially devised for stiff systems, the most popular being Gear's, 13 which allow a better time step regulation. Time steps of 10- 10 s are common when solving chemical kinetics problems with traditional schemes, and even methods specially developed for stiff systems usually do not allow time steps greater than 10- 6 s. Our steady-state problem clearly admits an equilibrium solution, which makes the use of time relaxation a possibility. However, we are solving for a large number of unknowns (typically about 500), and estimates of the characteristic species creation (or destruction) times and of the characteristic time scales of the fluid flow indi cate that we will be facing a stiff system. This implies a high price to pay in terms of CPU time when time relaxation is called, far greater than what a Newton method would require. Therefore, using this method throughout the whole computation seems inappropriate. However, it is a good way of bringing the initial guess, or the current approximate solution, back into the domain of convergence of the Newton method. Since the nonreac ting calculation preceeding the first call for time relaxation yields a flow field which is known to be close to the final solution for the reacting flow, we will consider the flow field to be constant when integrating. This means that it is not necessary to integrate the flow-field equations, which slightly reduces the CPU expense associated to time relaxation. The time integration routine used in the code is DDEBDF, developed by Shampine and Watts. 14 It is based on Gear's backward- differentiation formula,13 and allows solving a stiff system with accuracy up to order 8 with reasonable CPU expense. It is able to adjust the time step to minimize the CPU time expense while keep ing the accuracy of the computation within a limit specified by the user. A ma:<:i- mum number of ten calls to the time relaxation routine is allowed by REAC5 in its present version. The integration time is chosen by the user at the beginning of the calculation. In spite of its proven ability to integrate stiff PDE systems at a low cost, DDEBDF still requires high CPU expenses (a few days on a DEC VAXstation™ 3100) when applied to a typical CVD problem (say, HFCVD with 30 grid points and 25 species) . 32 The REA C5 code for CVD simulation It has been observed 17 that the introduction of a "scaling factor" S (typically 10- 3 or 10- 4 ) in front of the molar production rate Wk in the species conservation equation (3.1.12) makes it possible for DDEBDF to use a bigger time step. The effect of this scaling is to slow down chemistry in the system, making time derivatives less steep and thus inducing a lower local error on the solution. Of course, the use of S means that the total integration time should be multiplied by S to obtain the same result . In practice, it is not necessary to integrate this long, and a compromise can be found which allows integration over a reasonable time span (say 10- 3 s) and still brings the solution vector $ back in the domain of convergence of the Newton method after one or two integrations. 4.3.3 Modified Boundary Condition for the Mass Fractions Unlike the other boundary conditions, the mass flux balance on the substrate surface, Eq. (3.2.2), cannot easily be transformed into an unsteady condition. AI· though it is of the form Fi ( $) = 0 in the steady-state problem, it is not a conservation equation, and replacing the zero right -hand side by the time derivative of the wall mass fraction does not make sense physically. It is therefore necessary to modify this condition. Eq. (3.2.3) in section 3.2.3 can also be written: (4.3.1) Let us neglect thermal diffusion in the diffusion velocity. (Thermal diffusion is nonzero for the species H and H2 only, and since the output of the integration routine will be corrected by the steady-state Newton iteration which includes thermal diffusion, we can allow inaccuracy in the unsteady computation of the mass fractions). For the same reasons , let us assume: ~0 (OM) Oz W (4.3.2) Also notice that C k ,W is independent of time, since the wall temperature is fixed , that mk = Mk/N A and that, since we consider a unit volume, nk = NA[Xk], The REAC5 code for CVD simulation 33 and therefore, ank = NA a [XkJ at at = NAWk Under all these assumptions, Eq. (4. 3.1 ) becomes: After switching space and time derivatives in the second term of the left-hand side and using the finite-difference formula (4. 5.1 ), we can isolate the time derivative of the wall mass fraction: 1 'Yk M k -Ck Wk . +Dk,w ] aYk) - (aYk) -( at W - ------;,,---P u + D • .w [ I -:If 4 (W at W+l W W where (w = ZW+l - (w ZW· The right-hand side of this equation is known by the species conservation equation at the node VV + 1. The penalty to pay for use of this formula is a bigger Jacobian upper bandwidth, since the species conservation equation at the node W + 1 depends on variables at the nodes W to W + 2. 4.3.4 Modified Boundary Condition for Wall Axial Velocity As we saw in section 3.2.4, CVD induces a non-zero axial velocity Uw at the substrate surface ("bulk velocity"), whose expression is given by Eq. (3.2.4). When deriving the modified (unsteady) boundary condition for the wall mass fractions, we assumed the bulk velocity to be nearly independent of time (see Eq. (4.3.2)). This means that Eq. (3.2.4) should be replaced by the condition: However, changes in the mass fractions can lead to a change in the bulk velocity, which should be taken into account. In REAC5 , this is done by periodically exiting the time integration routine, updating UW with help of Eq. (3.2.4), then resuming integration. The "intermediate exit time" is chosen by the user at the beginning of the calculation, and satisfactory results have been obtained with an intermediateexit time equal to one tenth of t he total integration time. 34 The REAC5 code for CVD simulation 4 .4. Initial Guess for the Newton Method 4.4.1 Flow Field The majority of HFCVD cases involves flows at low Reynolds numbers (much smaller than 1). Therefore, it seems natural to take the creeping flow approximation as initial guess for the flow field (Goodwin and Gavillet 15 ) . In the case of a nonrotating substrate, and neglecting convective terms with respect to viscous terms in Eq. (3.1.13), (3.1.14), (3.1.15) and (3.1.16), we obtain: A V(z) = - - z(H - z) 2J.L T (z) = TL ~ + Tw (1 - ~) For a density profile p(z ) of general form, there is not necessarily an analyt ical solution for the axial velocity profile u(z). On the contrary, if the density gradient is approximated by the expression: '!:. dp ~ UL (TL - Tw) ~ pdz H TL H and if we define: T= TL -Tw 2TL we then find the approximate solution: as well as an approximate value for the radial pressure gradient parameter: OJiuLT H3 (4.4.1) Test simulations have shown that this initial guess lies within the domain of convergence of the Newton method for the first part of the computation (flow field calculation) even in cases where the creeping-flow assumption is not warranted. Therefore, the flow field profiles determined in this section are used as initial guess whatever the inlet conditions are. In case of a rotating substrate, the initial guess on the circumferential velocity is: The REA C5 code for C VD simulatIOn 4.4 .2 35 Pressure Gradient Parameter This approximated value of A is valid for a flow on a nonrotating substrate, which has to be "forced" by imposing a negative (in our coordinate system, defined in Fig. 3.1. 1) radial pressure gradient. In the case of a rotating substrate, however , it is possible for a flow to occur in the absence of a radial pressure gradient (due to the "pumping effect" of the surface, as in the von Kirmin 16 solution for a viscous stagnat ion flow on an infini te rotating plate in a semi-infinite domain) or even with an adverse (positive) gradient. The sign of the pressure gradient is also the sign of A. The adjustment technique for A described in section 4.1.1 does not allow a sign change. Therefore, in the case of a rotating substrate, it is vital to know beforehand whether the flow to be simulated is "forced" (A < 0), "ideal" (A = 0), or "starved" (A > 0). At a given substrate angular velocity, there is a one-t o-one relationship between the radial pressure gradient (and therefore A) and the quantity of fluid "pumped" by the rotating plate (and therefore the inlet axial velocity). A zero pressure gradient corresponds to an "ideal" inlet axial velocity, and it is clear that any inlet velocity larger (in modulus) than this ideal value will yield a forced flow, whereas any lower value (in modulus) will yield a starved flow. This fact is used in the code: if the specified substrate angular velocity is nonzero, the parameter A is set to zero and the ideal solution is obtained. The moduli of ideal and user-specified inlet axial velocities are then compared and the sign of A is chosen accordingly. There remains t he task of picking a reasonable initial value for IA I. Substrate angular velocities in CVD experiments are typically of t he order of 1000 rpm (as in Coltrin et alY for silicon CVD ), which renders the creeping-flow assumpt ion invalid . However , test simulations have shown that the modulus of the Ao defined in (4.4.1 ) can be used as initial value and still allow convergence. A slightly refined guess , used in the present version of REAC5 , is O. l lAol for a starved flow and - O.OOlIAol for a forced flow (these constants were determined empirically ). The REAC5 code for CVD simulation 4.4.3 36 Chemical Species Diamond CVD is often characterized by super-equilibrium concentrations of numerous species near the substrate surface. During the development of the code, it has been noticed that taking as initial guess a linear interpolation between substrate and inlet equilibrium mass fraction is a bad choice, since it is nearly always outside the domain of convergence of the Newton method. A more helpful solution is to define a wall residence time tre •. W and, st arting with the inlet composition, to integrate the rate equations at the temperature Tw over this time. The residence time tre..w can be a characteristic time scale of the flow . In the present version, an "average" inlet Peclet number is computed (that is, with an average mass diffusivity), and the residence time is chosen to be either H /U L or H2 / D depending on whether P e > 1 or not. Once the substrate mass fraction s are obtained, a linear interpolation is performed between inlet and substrate . 4.5. The Computational Mesh The choice of the computational mesh results from the classical dilemma between accuracy and cost of computation. On one hand, CVD diamond synthesis requires an accurate description of the situation near the substrate, especially as far as mass fractions and temperature gradients are concerned. This means that the spatial resolution near the substrate should be high. On t he other hand, the bigger the number of mesh points, the larger the system of equations (and therefore the Jacobian matrix) will be. This means that the computation will become more and more expensive. (For example, the number of Jacobian coefficients to be computed grows linearly with the number of points.) A uniform grid is unable to provide a good compromise. Fortunately, there are two ways to solve this problem: either operate on a fixed nonuniform grid, or use an adaptive, gradient-resolving gridding. The REAC5 code for CVD simulation 4.5.1 37 Fixed nonuniform grid The use of a fixed nonuniform grid, although it is not as "flexible" as adaptive gridding, is a simple and efficient way of concentrating points in a given region of the computational domain while keeping the total number of points, and therefore the CPU expense, within reasonable limits. The possibilities for grid generation are infinite. Keeping in mind that we want to concentrate points near the substrate surface and to be able to easily modify the grid, it seems natural to look for an analytical, parametric transformation which would change a uniform grid into one with higher node density near the origin z = O. Let us change the actual coordinate z into a new coordinate x by the transformation x = Xa (z) such that x H 1- e a where H is the reactor height and a is the transformation parameter. (Note that this transformation is not defined if a = 0, but that lima_o xa( z) = z) . This leads to a transformation of first and second partial derivatives In the equations (3 .1.12) to (3.1.17), according to the following rules: d dx d -=-dz dz dx a - e- a'71d = -a 1 - edx and ~ (F~) - (dx)2 [dF ~ + F!:....] + a2x2 F~ dz dz dz dx dx dx 2 dz dx 2 2 -a' ) 2 [dF d F d ] a -air F d = ( l_e-a e 71 dxdx + dx2 - H ( I_ e- a)e dx a The first and second derivatives of x with respect to z are fixed quantities, which are computed analytically at the beginning of the calculation. Consider now an uniform grid in the new coordinate x. If we take a > 0, it is the image by this transformation of a grid which is thicker near the surface z = O. The concentration of points varies according to the value of a (see Fig. 4.5.1 ). The choice of an uniform grid in the x-coordinate makes the programming easier. The The REAC5 code for CVD simulatIOn 38 1.0 0.8 0.6 :::r: "X 0.4 0.2 0.2 0.4 0.6 0.8 1 z / H Figure 4.5.1: Parametric coordinate transfonnation xa(z). Values of a are -2,-1,-0.5, 0.5, 1, 2. parameter a is fixed by the user at the beginning of the computation. Values of a between 0 and 2 have been used in the work presented here. 4.5.2 Adaptive gridding The coordinate transfonnation method has the inconvenience of concentrating grid points in a fixed way. Indeed, the choice of the transfonnation function allows to concentrate points near the substrate (a > 0) or near the inlet (a < 0). The problem is then to pick the value of a which best suits the case to be studied. For example, if the chemical mechanism is such that production of some species are important outside the region where points are concentrated, then the calculation will be less accurate and convergence will be more difficult. Therefore, the choice of the grid requires some prior knowledge of the problem to be solved. For the code to be flexible , and to adapt easily to a new problem, the choice of the grid should be problem-independent. Hence, it would be interesting to be able to adjust the The REA C5 code fOT CVD simulation 39 grid so as to concentrate points where needed and nowhere else. This can be done by using an adaptive grid. There are numerous criteria for grid adaptation. In our case, the key to accurate computation is a good modeling of the situation at the wall, which depends on a good calculation of the diffusion velocities Uk ,W. Since they are dependent on temperature and mass fraction gradients, the objective of grid adaptation should be to concentrate nodes in the regions of steepest variation of temperature and mass fractions. Therefore, it is natural to pick gradient resolution as the criterion for grid adaptation. We require the difference between any unknown G at two consecutive nodes i - I and i to be smaller than a sp ecified fraction {, of the maximum variation of this unknown over the computational interval, i.e. : maxC,i IGi - Gi-ll < {, I maxi G - mini GI This inequality is checked at each node i. If it is not satisfied, a new node IS inserted between i and i + 1, an initial guess is made at this node by linear interpolation, and the computation continues. In the code, a maximal grid size of 101 points cannot be exceeded. In our problem , flow-field variables and chemical variables (i.e., mass fractions ) vary over totally different ranges: a mass fraction can easily vary over ten orders of magnitude, whereas a flow-field variable does not. If the same tolerance {, is chosen for the mass fractions and flow-field variables, it will lead either to an inaccurate calculation of the flow field (if {, is chosen to be large, to accomodate mass fraction variations) or to an excessive accumulation of nodes in regions of steep mass fraction gradients (if {, is small, to accomodate flow-field calculation), which results in an unwanted additional CPU expense and possible round-off errors inducing uncontrolled oscillations of some solution components. A way around this problem is to take two different tolerances, one for the flow-field variables (usually {, fl ~ 15%) and one for the mass fractions ({,ch ~ 30%). Finite-differencing on a nonuniform grid is not as simple as on a uniform grid. Let us compute the approximate first and second derivatives of t he functi on f ( z) at an "inner" node i (i.e. not at the boundary of the computational domain). Consider The REA C5 code for CVD simulation 40 the points Zi_l, Zi, Zi+l , where I has the values Ii-I ,ii, };+I · Let us perform a parabolic interpolation of I. If I (z) ~ bo + biz + /),zz\ we find: 1 1 };-I bl = -yo 1 I. > 1 Ii+1 1 Zi-l zL z2 2' zi+! and /),z=~1 1 Zi Zi+l where T is the van der Monde determinant: 1 Zi-I Zi Zi+1 T = 1 zr-I II zl - (Zk - ZI) zl+1 '+I~k>I~'-1 The first derivative at the node i is then given by: 1 (:~) i ~ 2/),zz. + bl f (. - (.-1 f )'+ I + I" I" • ~i "'i + ~i-l ~i~i - l ~ I" (I" (i-I I" (. where (. = Zi+1 - Zi· Similarly: (~{). ~ (1.:::: =~~;i-t 2 2 2 ~ (. ((. + (.-1) };+I - ((i-I Ii + ((i + (i-I) (.-1 Ii-I where the intermediate derivatives (index i ± 1/2) are computed using a linear interpolation of the type: dI) Ii+1 - Ii ( dz .+! - Z'+1 - Z. At the first node (wall ), only first derivatives are needed in the computation of the diffusion velocities. They are approximated by forward linear interpolation: dI) _ IW+I - Iw ( dz W ZW+I-ZW (4.5.1) First derivatives at the last node are approximated by a backward linear interpolation between the nodes N m - 1 and N m , analogous to Eq. (4.5.1 ). The REAC5 code for CVD simulation 41 4.6. Problems Specific to the Last Node We have seen that the boundary conditions do not allow a direct calculation of the parameter A, and that therefore, it has to be iteratively adjusted. Imposing a value of A determines the pressure gradient, and has a direct infiuence on the inlet axial velocity UL . Therefore, the condition on the inlet axial velocity is released and the value UL becomes part of the calculation. This means that the axial momentum equation has to be solved at the inlet node too. Since we are at the boundary of the computational domain , we cannot use centered differencing. A full backwarddifference scheme on this second-order equation would involve variables from the last two nodes preceding the inlet node. Since the Jacobian matrix is stored in banded form, this would impose a larger lower bandwidth for the whole Jacobian matrix and the storage of useless zeroes at all nodes but the last. A way to avoid this difficulty is to create a fictitious node somewhere between the last two, and to solve the axial momentum equation at this point using centered differences. A value Zin t of the axial coordinate (in the present version of the code, Zint = O.lZNm_1 + O.9ZNm ) is picked in [ZNm-l, zNml . A parabolic fit lS of viscosity, density, axial and radial velocity is performed over the last three nodes (Nm - 2 to N m ), and is used to evaluate these variables at Zint. The derivatives at the intermediate node are computed using the formulae of section 4.5.2, the only difference being that Zi, fi replace Zi _ l, !i- I and Zint, fint replace Zi, fi in all expressions. 42 Chapter 5. Results 5.1. Validation of the Code 5.1.1 Flow Field Calculation The nonreacting part of the calculation has been tested against a reacting computation performed by Coltrin et al. ll for the study of silicon CVD , with a code similar in principle to REAC5 , except for the fact that the axial momentum equation is not solved (the axial velocity is determined from the continuity equation) and that the static pressure is assumed constant in the z-direction. Silicon CVD is similar to diamond CVD: a gaseous feedstream of hydrogen and silane (Si H4 ) with more than 99% hydrogen is injected at a low pressure and impinges on a silicon substrate. However, injection occurs at or near room temperature, and it is the substrate that is heated. Cases involving both rotating and nonrotating substrates have been investigated. Unfortunately, data in CHEMKIN-II form about silicon compounds SiHj (i=l ,2·, i=0,4) , namely coefficients of the polynomial fits of the thermodynamic properties, were not available, which prevented calculation of species profiles as well. The results presented by Coltrin et al. are for a heated-plate-type flow: a mIxture of silane and hydrogen (99%H2/1 %SiH 4 ) diluted in an overwhelmingly present carrier gas (argon) is injected at 300 K and atmospheric pressure, and impinges on a substrate heated to 1000 K. The reactor height is 1.1 cm. Four cases (three rotating and one nonrotating) have been simulated. Whereas nonrotating reactors require a negative radial pressure gradient for a flow to occur, rotating reactors can be "ideal" (zero pressure gradient), "forced" (a favorable pressure gradient is added, which results in a higher quantity of pumped gas and therefore an inlet velocity of higher magnitude), or "starved" (the flow is hindered by an adverse pressure gradient, which has opposite effects) . All three cases have been investigated. Results 43 With a substrate angular velocity of 1000 rpm and a pure argon feedstream , the ideal inlet axial velocity predicted by REAC5 is approximately -8.59 cm/s. The value computed by Coltrin!! with his reacting argon/hydrogen/silane feed stream is -9.01 cm/s. This, combined with the close resemblance of velocity and temperature profiles, leads to the confirmation of two hypotheses. First , the differences in the transport properties of the actual feedstream and the pure carrier gas are small (as assumed, for example, by Evans and Greif!2 ) and the flow field can be predicted fairly accurately when considering the carrier gas alone. Second, the flow fields in reacting and nonreacting cases differ very little. The starved case was chosen with an inlet velocity of -2.3 cm/s, and the forced case with -23 cm/s. Fig. 5.1.1 , Fig. 5.1.2 and Fig. 5.1.3 show good agreement of the REAC5 flow field profiles with Coltrin's. The nonrotating simulation was performed with an inlet velocity of -115 cm/s. This high value induces boundary layer-type radial velocity and temperature profiles, which have been accurately predicted by REAC5 (cf. Fig. 5.1.4). A comparison of four cases involving low and high inlet velocities, rotating and fixed substrates, positive and negative radial pressure gradient, shows no major discrepancies between the results of REAC5 and those of Coltrin et al. 11 The slight differences in the profiles can be explained by the difference in feedstream composition, the nonreacting nature of the REAC5 simulation compared to the reacting Coltrin simulation, and by the difference in solution principle between the codes (5 flow field equations and variable pressure versus 4 equations and fixed pressure). Therefore, we can consider that the flow field calculation is validated. 5.1.2 Species Profiles Calculation It has been previously explained that a simulation of Coltrin's silicon CVD experiment including reactions could not be performed due to the lack of data concerning silicon-hydrogen compounds. Also, no report has been made of the application to diamond CVD of any of the known simulation codes. Therefore , it was not possible to perform a priori validation of the species profile calculations. However, comparison of simulation results with diamond CVD experimental data, Results 44 Figure 5.1.1: Flow field for Coltrin's "starved" case. (Dashed curve is Coltrin 's calculation.) as will be shown, yielded good agreement, taking into account the simplicity of the surface chemistry model and the three-dimensional nat ure of the flow in the actual reactor. 45 Results 0 -2 15 -4 < ., "E 10 -6 -., "- ~ 0 \ ::> -8 \ / 5 ~ ~ -10 -12 ~ 0 0.2 0.4 0.6 0.8 1 80 0.2 0.4 0.6 0.8 1 800 -; 60 A "...... 3 0 1000 100 ., 0 600 40 20 0 .. 400 0 0.2 0.4 0.6 0.8 1 X ( cm) 0 " 0.2 0.4 0.6 0.8 1 X (c m) Figure 5.1.2: Flow field for Coltrin's "ideal" case. (Dashed curve is Coltrin's calculation.) 5.2. Hot Filament CVD Simulation T he code has been applied to simulation of two HFCVD experiments for which a detailed description of the operating conditions are available. Predictions of species molecular density and diamond deposition rates have been compared with experimental data. 46 Results 30 0 25 -5 20 < ~ -10 en "E u 15 ~ " " -20 0 10 '< "" 0 .2 0.4 0.6 O.B "en ~ ~ -15 :;, ~ ~ ~ ~ 5 1 0 " 0 0.2 0.4 0.6 O.B 1 1000 100 80 BOO -; 60 en A "- 600 40 :3 20 0 ~ 400 0 0.2 0.4 0.6 O.B 1 X (em) 0 0.2 0.4 0.6 O.B 1 x (e m) Figure 5.1.3: Flow field for Coltrin's "forced " case. (Dashed curve is Coltrin's calculation.) 5.2.1 Description of First Experime nt In the experiment of Harris et al.,19 a feedstream in the proportions of 0.29 sl/min CH 4 for 100 sl/min of a gas consisting of 1000 ppm of Ne in H2 is injected at room temperature past a tungsten carbide filament at 2600 K, and impinges on a silicon substrate. The distance between filament and substrate can be adjusted, 47 Results 800 600 400 o 0 .2 0.4 0 .6 0 .8 1 X (em) Figure 5.1.4: Flow field for Coltrin's nonrotating case. (Dashed curve is Coltrin's calculation. ) and results are given for values of H up to 4 cm. The static pressure at the inlet is 20 Torr. According to Harris, optical pyrometer measurements show that the gas temperature as it leaves the filament is approximately 2000 K, and that the substrate temperature is approximately 1000 K. Moreover, measured concentrations of the principal species (methane, acetylene, methyl radical), compared with results of Results 48 numerical integration of the rate equations (3.2.1) at the constant temperature 2000 K, suggest a "residence time" of t he gas near the filament of about 0.2 s, which is not sufficient to reach full chemical equilibrium. Simulations of this experiment were performed with a fixed (non-adaptive) grid , for the reactor heights 0.5, 1, 2 and 4 cm. The inlet velocity was -0.5 cm/s except for the case H = 0.5 cm where the two values -0 .5 and -1 cm/s were investigated. The grid was nonadaptive and uniform (transformation parameter a = 0), and contained 31 points in all cases except H = 0.5 cm and UL = -1 cm/s, where it contained only 26 points. As soon as adaptive gridding became available, comparison tests were performed, which showed only minor differences in the calculated species profiles and diamond growth rates, but a significant improvement in convergence time (5 hr 10 min vs. 5 hr 50 min , i.e. , 10.3%, on a DEC VAXstation™ 3100). 5 .2.2 Description of Second Experiment The experiments conducted by Butler et al. 22 have also been investigated. The feedstream composition is 96 to 99.9% H2 and 4 to 0.1 % CH 4 , injected past a tungsten filament at 2500 K. The substrate is also silicon, and its temperature is 1200 K. The reactor height is 2 cm, and the inlet static pressure is 20 Torr. However, we assumed the same gas residence time and inlet temperature as III Harris 's experiment, since both filament temperature and feedstream composition were similar. This means that we considered an inlet temperature TL = 2000 K and an inlet residence time of 0.2 s. Results include measurements of species number densities 0.5 cm away from the substrate surface for a 1% CH 4 feedstream , obtained by infrared diode laser absorption (IR) , resonance-enhanced ionization (REMPI) or laser-induced fluorescence spectroscopy. Data is available for the species CH 4, C 2H 2, C2H 4 and CH 3, and an upper bound has been determined for the species CH2, C2H6 , C3H4, C3H6 (methylacetylene and cyclopropane), C3HS and C 4H2. Simulations were performed with a nonadaptive uniform grid featuring 21 (case H = 2 cm, 0.5 % CH 4) or 31 points (other cases), for the reactor heights 0.5, 1.5 and 2 cm, and the methane percentages 0.5, 1, 2 and 4 %. 49 Results 5.2.3 Chemical Mechanism for Simulation The mechanism describing the chemical process inside the reactor has been chosen to include all the C 1 and C2 hydrocarbons except C and C2 , all the C3 except C 3Hg and some of the C 4. C 1 and C 2 chemistry is based on the mechanism used by Harris 23 in previous CVD experiments. However, our mechanism includes pressure-dependent rate constants for the unimolecular reactions describing pyrolysis of CH 4, C2H6, C2Hs and C2H3. The limiting high-pressure rate constants were determined by Warnatz,24 and we assume a Lindemann fall-off function. Chemistry for the C 3 and C 4 hydrocarbons considered is taken from Harris et al ..2S The mechanism file , in a form suitable for use in the CHEMKIN-II interpreter,S can be found in Appendix B. It includes a total of 2 elements, 25 species and 56 reactions. The same mechanism was used in the simulation of both experiments. Due to lack of thermodynamic data in CHEMKIN-II form for neon , the actual feedstream used by Harris has been replaced by a mixture of hydrogen and methane in the same proportions. The inlet composition was determined by integrating the rate equations during 0.2 s at the constant temperature 2000 K, using the mechanism previously mentioned. 5.2.4 Results The results of these two experiments and of their numerical simulations have been discussed by Goodwin and Gavillet,15 and we will briefly recall the main conclusions of this discussion . • First Experiment (Harris et al.) Growth rates between 0.5 and 1 pm/hr have been reported. The other available experimental results are mole fractions of the major species (methane, acetylene) at the substrate surface. The reactor geometry is unfortunately not close to the ideal one-dimensional geometry of the code: the reactor height is not much smaller than the subst rate diameter, the filament is linear instead of covering the whole injection area, and does not span the entire width of the injection section. Therefore, we should expect noticeable discrepancies between experimental data and simulation results. Methane Results 50 and acetylene mole fraction measurements at the substrate surface show that the methane mole fraction approaches the inlet (unreacted) mole fraction as the distance between filament and substrate is increased. This strongly hints to diffusion of cold, unreacted gas into the warm region below the filament , a factor that cannot be included in our one-dimensional model. The flow field in the operating conditions of Harris et al. is very close to the corresponding initial (creeping flow) guess. This is hardly surprising, since the height-based Reynolds number is very small (e.g., with H = 0.5 cm and UL = -1 cm/s, Re = 5.33 x 10- 4). Fig. 5.2.1 shows a comparison between initial guess and actual result. The possible role of several hydrocarbons in the diamond CVD process has been investigated. To get an upper bound of the contribution of a species to film growth, its reaction probability I has been set to 1 and all others to O. Methyl and acetylene, the proposed zo .Z! precursors of diamond formation, have been treated differently: their reaction probability has been adjusted so that each of these species, taken separately, yields a growth rate comparable to all others in the conditions of Harris 's experiment. The results are listed in Table 5.2.1. It should be kept in mind that a species yielding a high growth rate does not necessarily contribute significantly to diamond growth; it simply cannot be ruled out on the basis of gas-phase chemistry arguments. On the other hand, a species yielding a low growth rate is almost certainly not a precursor of diamond growth. First of all, these results confirm the possible role of the methyl radical and acetylene in the deposition process: with reaction probabilities at least three orders of magnitude below those of other species, they cause similar growth rates. They are the only species which could realistically account for the 0.5 to 1.0 /l-m/ hr growth rates obtained by Harris et al.. More surprising is the fact that the heavier species C2 H6, C3HZ, C3H3 and C3H4 , and, to a lesser extent, the species C3Hs and CHz, could contribute to diamond growth (although their contributions with I = 1, in regard of the reported growth rates, make it unlikely). However, there is no production of C 3R2 near the substrate, Results 51 1.5 2000 i' /' ,..... / 1.0 V /' /' / :-.... / :---.. /' / ---> en / '/ 1 u ~ ~ / / r E ~~ // T /' 0.5 1800 / ~ 1600 / / 0 I- / ::J 1400 / / .... ~ / -0.5 -~ / ~ ~ u / / ~ ~ / / ~ ~ .... 1200 / -1.0 ~ 0 0.1 0.2 0.3 0.4 1000 0.5 z (em) Figure 5.2.1: Comparison between creeping flow approximation (dashed line) and calculation result (solid line) for H = 0.5 cm and IS tiL = -1 cm/s. (From Goodwin et al. ) and its net production in the whole reactor is negative. 15 This makes it unlikely for C3H2 to be a contributor to growth. On the contrary, the other species of Table 5.2.1 all have net positive productions. Even if a surface reaction probability of 1 is somewhat unlikely, this computation shows that these heavier hydrocarbons , and especially C3H4 , cannot be ruled out as minor contributors. 52 Results Table 5.2.1: Calculated film growth rates. Species Growth rate (/lm/hr) C H2 1 0.075 CH 3 10- 3 0.14 C2 H 2 7 x 10-5 0.1 C 2H 5 1 0.022 C2 H 6 1 0.117 C3 H 2 1 0.235 C3 H 3 1 0.156 C3 H 4 1 0.232 C3 H5 1 0.093 • Second Experiment ( Bu.tler et al.) As in the experiment of Harris et al.,19 the reactor used by Butler et al. 22 is not fully axisymmetric, and a non-negligible fraction of the feedstream is not heated enough by t he filament to react . Therefore, the same restrictions on the accuracy of concentration measurements apply. The observations made on the flow field in the previous section apply here too, for the same reasons. The number densities reported by Butler et al. agree reasonably well with the values predicted by REAGS (see Table 5.2.2). The validity of the experimentally determined upper bounds is confirmed. The biggest discrepancies concern methane and acetylene concentrations. It should be noticed that Butler's values are pathlength-averaged over the whole chamber, which is 25 cm wide. The fact that the filament only covers a small fraction of the injection area causes a nonnegligible fraction of the injected gas to remain at or near room temperature, which prevents methane decomposition and formation of other hydrocarbons. This means that the averaging process inherent to the measurements takes into account both hot and cold regions, and that therefore, the experimentally measured methane Results 53 Table 5.2 .2: Measured and calculated species number densities. Species Measured (Ref. 22) Calculated (cm- 3) (cm - 3) 2 ± 1.5 x 10 12 8 ± 3 x 1014 1.6 X 1013 1.1 X 10 14 2 ± 1 x 10 14 6 ± 2 x 10 12 6.0 X 10 14 2.0 X 10 12 9.4 X 10 10 C 3H 6 < 8 X 10 13 < 1.4 X 10 14 C 4H2 < 3 X 10 13 2.1 X 1011 CH 3 CH 4 C2 H 2 C2H4 C 3H4 2.2 X 109 concentrations should be higher (and the methyl concentrations lower) than the values predicted by the code. The results confirm this hypothesis. It should also be noticed that the simulations were performed with all surface reaction probabilities set to zero, but it has been observed in the simulations of the Harris experiment l9 that concentrations a distance H /4 away from the substrate do not depend strongly on the reaction probabilities. 5.3. Preliminary Plasma Torch CVD Simulation Provided that the temperature does not exceed 5000 K (upper bound of the domain of validity of the polynomial fits performed by CHEMKIN-II on thermodynamic data), it is possible to simulate plasma torch CVD with the REAC5 code . A preliminary simulation of an experiment by Cappelli et al. 26 at Stanford University has been performed. Results 5.3.1 54 Experiment The plasma torch of Cappelli et al. transforms a 157/ 16/1 mixture of ar- gon,hydrogen and methane (molar proportions) into a plasma by exposure to a radio frequency (RF) field. The freestream Boltzmann temperature has been estimated to be 4900 K, and will be considered to be the inlet temperature in the simulation. A water-cooled substrate is placed 8 em away from the torch nozzle. However, according to Owano,27 the gas flow becomes uniform only 2 em outside the torch. Therefore, we will retain 6 em as reactor height, the other parameters being unchanged. The substrate temperature has been estimated to be 1400 K in these conditions. With a total gas flow rate of about 120 sl/min (at 1 atm and 290 K ) and an internal torch diameter of 7 em, the axial velocity at the nozzle, which will be considered the inlet axial velocity in the simulation, is approximately -913 cm/s. The simulation was performed with an adaptive grid, initially featuring 26 unevenly distributed points (t ransformation parameter a = 2). The final grid featured 30 points, 2 of which were added during the nonreacting part of the calculation to resolve the important gradient of V near the wall. 5.3.2 Chemical Mechanism for Simulation In our preliminary simulation, the operating conditions, including feedstream composition, were identical to Cappelli's. However, the only chemical reaction taken into account was hydrogen recombination. The purpose of this calculation was to confirm experimental observations of super-equilibrium concentrations of atomic hydrogen near the substrate. The mechanism file contained 3 elements, 4 species and 1 reaction. 5.3.3 Results As expected, the flow field is of boundary-layer type: both radial velocity and temperature profiles show the existence of a velocity and of a thermal boundary layer. P eciet numbers are greater than 50, suggest ing a convection-dominated flow. The existence of a thermal boundary layer, along with the chemical equilibrium Results 55 assumption at inlet, explains the absence of chemical production over a large frac tion of the reactor height where the temperature is constant. I I 1000 r- I I I - '" '- ('") •• E '" o ... '- <> L o ... E '" - .c - ... ::> o ., -1000 - - -2000 - - '"Eo • ·• <D 0 r I I I I I 0 1 2 3 4 5 6 x (em) Figure 5.3.1: Molar production rates, considering only hydrogen recombination. The influence of diffusion can be seen in this case. Molar production rate profiles (Fig. 5.3.1) show destruction , then production of molecular hydrogen as the gas enters the thennal boundary layer. This is confirmed by a comparison of the temperature profile with the profile calculated for a nonreacting mixture in the (fixed ) inlet composition (Fig. 5.3.2): the reacting profile undershoots, then overshoots the nonreacting profile as the substrate is approached, due to the heat release associated with hydrogen recombination. Yet, mass fraction profiles for atomic and molecular hydrogen vary monotonically between inlet and substrate (Fig. 5.3.3), instead of showing extrema. This can been explained by diffusion of atomic hydrogen towards the substrate. Results 56 4000 ~ ::> co rt 2000 o 2 3 4 5 6 X (em) Figure 5.3.2: Comparison of temperature profiles for nonreacting flow and flow with hydrogen recombination. It could be interesting to check whether the reaction is locally equilibrated or not. The variable to be considered for this purpose is the partial equilibrium ratio: where Kc is the equilibrium constant of the two-way, three-body reaction: 2H + M ;=' H2 + M in terms of molar concentrations (M is any other molecule). If E > 1, molecular hydrogen is in super-equilibrium concentration, whereas atomic hydrogen is if E < 1. Fig. 5.3.4 shows the evolution of E as the gas approaches the substrate. Starting from equilibrium at the inlet, the gas remains in chemical equilibrium until well into the thermal boundary layer (Fig. 5.3.2 shows a boundary layer thickness of about 2 cm, whereas E differs from 1 only for z < 1.2 cm). Then, E increases to reach a peak value of approx. 1.15 at z c:: 0.55 cm, before dropping sharply below 1 ( E = 1 at z c:: 0.25 cm) to reach 1.59 x 10- 8 at the wall. This means that molecular hydrogen 57 R esults JI I I I I A. - - .- " ~ ::J rt CII. - - -5~~~~I~~.~I~~~I~~~I~~~I~w o 1 2 3 4 5 6 x (em) Figure 5 .3 .3: Mass fractions , considering only hydrogen recombi- nation. is in super-equilibrium concentration between z ~ 0.55 and z ~ 0.25 cm, whereas it is atomic hydrogen which is in super-equilibrium concentration below z ~ 0.25 cm. The presence of an excess of molecular hydrogen in the central part of the thermal boundary layer is consistent with its destruction in this region (an attempt by the system to bring itself back to equilibrium), and the origin of this excess could be a quicker diffusion of atomic hydrogen towards the substrate surface. The higher thermal diffusion ratio of atomic hydrogen would tend to confirm this hypothesis. 58 Results I I I I I 1.00 - - 0.75 .,<D W . ~ ::> .,'" L rt "::> 0.50 - 0.25 - .0 " I I I I I . 0.00 ~....L.....l-'-:~....L.....l-'-:~....L.....l-'-:~....L.....l-'-L........L.....l-'-L......L.....l~ o 1 2 3 456 x (em) Figure 5.3.4: Partial equilibrium ratio for hydrogen recombination. Chapter 6. Further Improvements The REACS code seems to work reasonably well in its present version. However , a certain number of improvements can be suggested. First of all, the current surface chemistry model, with its surface reaction probabilities, is a very crude and empirical one, and should be replaced by a more realistic model. The 25-species, 56-reaction mechanism currently used to describe methane pyrolysis could be refined. For example, when operating at high temperatures (e.g. plasma CV D), it is likely that dissociation is going to occur on a large scale, and therefore the mechanism should include atomic carbon as a species. Unfortunately, data on reactions involving atomic carbon are sparse. Convergence problems have arisen during plasma CVD simulations, which feature boundary-layer-type temperature and velocity profiles and high inlet temperat ures. The origin of these difficulties could as well be the nature of the problem itself as an initial guess on the species mass fractions too far from the actual profiles. At the present time, time relaxat ion is more to be avoided than sought due to its very high CP U cost . Any technique allowing a reduction of the CP U t ime required by time relaxation should be implemented. There is little hope for significant improvement in this area, due to the stiffness of the system and the lack of availability of an integration scheme significantly more performant than Gear's . REAC S has not yet been tested in combustion torch CVD cases. Due to the far greater number of species and reactions involved in a combustion mechanism and to the relatively high velocities at the torch nozzle, convergence difficulties and a high CPU expense may be encountered. 60 Chapt er 1. Conclusion A simplified one-dimensional numerical model calculating the reacting axisymmetric flow in various types of CVD reactors has been developed and applied to the simulation of hot-filament diamond CVD experiments. Preliminary simulations of RF plasma torch CVD have been performed, and confirm super-equilibrium concentrations of atomic hydrogen near the substrate surface. Although the current model has some shortcomings, especially in modeling the deposition process itself, the agreement with experimental data can be considered satisfactory. The code has provided insight on the nature of the flow field in the reactor , and has allowed us to assess the role of various hydrocarbons in the diamond synthesis process. 61 References 1. P.K. Bachmann & R. Messler, C & EN Special Report , May 15, 1989. 2. K.E. Spear, J.Am. Ceram. Soc. 72, [2],171-191 (1989). 3. J.C. Angus & C.C. Hayman, Science ill, 913 (1988). 4. L.M. Hanssen, W.A. Carrington, J .E. Butler & K.A. Snail, Mat. Lett. 1, 78,289-292 (1988). 5. R.J. Kee, F .M. Rupley, J.A. Miller, CHEMKIN-II: A Fortran Chemical Kin etics P ackage for the Analysis of Gas-Phase Chemical Kinetics, Sandia National Laboratories, report # SA D89-8009 (1989). 6. T. Eversole, Synthesis of Diamond, U.S. patent # 3030188, April 17, 1962. 7. B.V. Derjaguin & D.B. Fedoseev, Scientific American 233, [5],102-109 (1975). 8. S. Matsumoto, Y. Sato, M. T sutsumi & N. Setaka, J. Mater. Sci. 17, 3106-3112 (1982) . 9. S. Matsumoto, Diamond and Diamondlike Materials Synthesis, 95-98, edit ed by G.H . Johnson, A.R. Badzian & M.W. Geis, Materials Research Society, Pit tsburgh , PA (1988). 10. J.F. Grcar, R.J. Kee, M.D. Smooke & J .A. Miller, Twenty-First Symposium (International) on Combustion, 1773, The Combustion Institute, Pittsburgh, PA (1986) . 11. M.E. Coltrin, R.J. Kee, J. Electrochem. Soc. 136,819 (1989 ). 12. G. Evans & R. Greif, J. Heat Transf. 109, 928 (1987). 13. C.W. Gear, Numerical Initial Value Problems for Ordinary Differential Equations, Prentice- Hall (1971 ). 14. L.F . Shampine & H.A. Watts, DEPAC - Design of a User-Oriented Package of ODE Solvers, Sandia National Laboratories, # SAND79-2374 (1979). 15. D .G. Goodwin & G.G. Gavillet, Numerical Modeling of the Filament-Assisted Diamond Growth Environment, to appear in J . of Appl. Phys. (1990). 16. T. von Karman, Zeit. f. Angew. Math. u. Mech.l, 70-97 (1921 ), §7. 62 17. R.J. Kee, private communication (1990 ). 18. D.G. Goodwin, private communication (1990). 19. S.J. Harris, A.M. Weiner & T.A. Perry, Appl. Phys. Lett.~ , 1605-07 (1988). 20. M. Frenklach & K.E. Spear, J. Mat. Res . .3., 133-140 ( 1988). 21. M. Tsuda, M. Nakajima & S. Oikawa, J. Am. Cham. Soc. 108, 5780 (1986). 22. J .E. Butler & F .G. Celii, Proc. 1st Conf on Diam. & Diam.-Like Films, 317324, J.P. Dismukes, ed., The Electrochemical Society, Pennington, NJ (1989). 23. S.J. Harris, J. Appl. Phys. 65, 3044-3048 (1989). 24. J. Warnatz, Combustion Chemistry, 197-360, edited by W.C. Gardiner, Jr. , Springer-Verlag, New York ( 1984). 25 . S.J. Harris, A.M. Weiner & R.J. Blint, Combustion and Flame 72, 91-109 ( 1988) . 26. T.G. Owano, C.H. Kruger & M.A. Cappelli, Diamond Synthesis in a 50 kW Inductively Coupled Atmospheric Pressure Plasma Torch, to appear in Journal of Materials Research (1990 ). 27. T.G. Owano, private communication (1990) . 28. D.G. Goodwin, private communication (1990). 29. S.A . Self, Notes on Boundary Conditions for Diffusion to Reacting Surfaces, 1988 (unpublished). 63 Appendix A. M ass Flux Balance at the Substrate Surface To find the boundary condition on the chemical species at the wall, it is necessary to perfonn a detailed mass flux balance, including surface deposition , at the substrate surface. The following derivation was performed by Goodwin,28 based on notes by Self. 29 Let q and 'D~ be, respectively, the creation and destruction rates of the species k due to surface reactions (in the eGS system, their units would be molecules/cm 2 I s). Let wk+ and I¥k- be, respectively, the molecular fluxes of the species k leaving and arriving at the surface (i.e. with z-component of their molecular velocity respectively positive and negative, in the coordinate system of Fig. 3.1.1 ). Let us assume that both fluxes are Maxwellian. Then: where nk+ and nk- are, respectively, the number densities of molecules k with axial molecular velocity positive or negative (at equilibrium, nk+= nk- = nU2 ). The net fl ux leaving the surface is: W% = I¥k+ + I¥k= q -'D% ( A.l ) (note that I¥k- < 0) Assume that the destruction rate is represented by a surface reaction probability -n' .T.'<-'1=-,1""'1 The molecules which are not destroyed by surface reactions have to be reflected , and therefore: 64 Hence: Dividing by Ck/2 yields the number density n k+. Also, from Eq. (A.1 ), we obtain: which yields n k- . - ,+ + n ,- , or.. Bu t n , k - nk k (A.2 ) Now let : Combining Eq. (A.2 ) and the definition of the diffusion velocity (3.1.7) yields: Ck,W ( 'Yk) ] ['Yk 4- + 1 - 2" "W Yk ,w Mk kBTw C. MPw k _ (1- ~) D~w [Mk0 k (d (~:T))w - (dC:Yk ) ) J 65 Appendix B. Chemical Mechanism File The following is the chemical mechanism file, in CHEMKIN-II fonn, used in the simulations of hot-filament CVD experiments. Each reaction line contains, in that order , t he constants A, (3 and E necessary to the computation of the equilibrium constant ]{p: ELEMENTS HC END SPECIES CH CH2 CH3 CH4 C2H C2H2 C2H3 C2H4 C2H5 C2H6 C3H2 C3H3 C3H4 C3H5 C3H6 I*C3H7 N*C3H7 C4H C4H2 C4H3 C4H4 C4H6 C4H8 H H2 END REACTIONS H + H + M = H2 + M 9.7E16 -0.6 0.0 !Al HARRIS CH3 + H (+M ) = CH4 (+M) 6.0E16 -1.0 0.0 !WARNATZ LOW /8.0E26 -3.0 0.0/ !A2 HARRIS CH4 + H = CH3 + H2 2.2E4 3.0 8800.0 !A3 HARRIS CH4 + CH2 = CH3 + CH3 1.0E13 0.0 0.0 !A4 HARRIS CH3 + CH3 = C2H5 + H 8.0E14 0.0 26500.0 !A5 HARRIS CH3 + H = CH2 + H2 7.2E14 0.0 15100 . 0 !A6 HARRIS CH3 + CH2 = C2H4 + H 2.0E13 0.0 0.0 !A7 HARRIS CH3 + M = CH2 + H + M 1.0E16 0.0 90600.0 !A8 HARRIS CH2 + H = CH + H2 4 .0E13 0.0 0.0 !A9 HARRIS CH + CH4 = C2H4 + H 6 . 0E13 0.0 0.0 !Al 0 HARRIS CH + CH3 = C2H3 + H 3 . 0E13 0.0 0.0 !All HARRIS C2H6 + H = C2H5 + H2 5.4E2 3.5 5200 . 0 !A12 HARRIS 66 C2H6 + CH3 = C2H5 + CH4 C2H6 (+M) LOW /1. OE19 C2H5 (+M) LOW = CH3 + CH3 ( +M ) 0 .0 0.0 4.0 8300.0 !A13 HARRIS 2.4E16 0.0 87416.0 !WARNATZ 68070.0/ = C2H4 + H (+M) /1. OE17 5.5E-l 2.0E13 !A14 HARRIS 0.0 39648.0 !WARNATZ 31049 . 0/ !A15 HARRIS C2H5 + CH3 = C2H4 + CH4 7.9Ell 0.0 0.0 !A16 HARRIS C2H4 + M = C2H2 + H2 + M 2.6E17 0.0 79400 .0 !A17 HARRIS C2H4 + M = C2H3 + H + M 2.6E17 0.0 96600.0 !A18 HARRIS C2H4 + H = C2H3 + H2 1 .5E14 0.0 10200.0 !A19 HARRIS C2H4 + CH3 = CH4 + C2H3 4.2El1 0.0 11200.0 !A20 HARRIS C2H3 + H = C2H2 + H2 2.0E13 0.0 0.0 !A21 HARRIS 1. 6E14 0.0 37976.0 !WARNATZ C2H3 ( +M) LOW = C2H2 + H ( +M) /3.0E15 0.0 32000.0/ !A22 HARRIS C2H3 + CH3 = C2H2 + CH4 7.9E11 0.0 0.0 C2H2 + M = C2H + H + M 4.0E16 0.0 107000 .0 !A24 HARRIS C2H2 + H = C2H + H2 6 . 0E13 0.0 23700.0 !A25 HARRIS I*C3H7 = C3H5 + H2 1.0E13 0.0 30000.0 !HOl N*C3H7 = C3H6 + H 1.0E14 0.0 37280.0 !H02 I*C3H7 = C3H6 + H 2.0E14 0.0 38720.0 !H03 C3H6 + M = C2H3 + CH3 + M 1.0E18 0.0 74000.0 !H06 C3H6 = C3H5 + H 1.0E13 0.0 78000.0 !H07 C3H6 + H = C3H5 + H2 5.01E13 0.0 3500.0 !H08 C3H5 = C3H4 + H 3 . 98E13 0.0 70000.0 !Hll C3H5 + CH3 = C4H8 1.34E13 0.0 0.0 !H12 M C3H5 + H = C3H4 + H2 1 .0E13 0.0 0.0 !H14 C3H5 + C2H4 = C3H6 + C2H3 3 .2E12 0.0 18000.0 !H16 C3H4 + H = CH3 + C2H2 2.0E13 0.0 2411. 0 !H18 C3H3 + H = C3H4 2.0E13 0.0 0.0 !H21 C3H2 + H = C3H3 6.0E12 0.0 0.0 !H24 C2H2 + CH2 = C3H3 + H 2.0E12 0.0 0.0 !H27 C2H5 + C2H3 = C4H8 8.9E12 0.0 0.0 !I07M !A23 HARRIS 67 C4H6 + H = C2H4 + C2H3 5.0Ell 0.0 0.0 !J07 C2H3 + C2H3 = C4H6 8.9E12 0.0 0.0 !J12 C4H6 + H = C4H4 + H2 + H 1.0E14 0.0 14500 .0 !KOl C4H6 + C2H3 = C2H4 + C4H4 + H 6.3E13 C4H4 + M = C2H2 + C2H2 + M 1.7E93 C2H2 + C2H = C4H3 0.0 -20.5 14500 . 0 !K02 139000 .0 !K07 3 . 9E7 1 . 66 -1600 .0 !K08M 1.6E13 0.0 25060.0 !K09 C4H4 + M = C4H2 + H2 + M 5.6E92 -20 . 5 139000.0 ! Kll C4H4 +H = C4H3 + H2 3 . 89E14 0 .0 14500.0 !K04 + K12 C4H4 + C2H = C2H2 + C4H3 8.0E13 0.0 0.0 !K05 + K13 C4H4 + C2H3 = C2H4 + C4H3 1. OE12 0.0 16300.0 !K06 + K14 C4H3 + H = C4H2 + H2 2.0E13 0.0 o. 0 ! LO 1 + L08 C4H3 + M = C4H2 + H + M 8.5E15 0.0 60000 . 0 !L02 C4H2 + M = C4H + H + M 3 . 5E17 0.0 80000 . 0 !L03 C4H + H2 = H + C4H2 4 . 0E13 0.0 0.0 !L06 4 . 0E13 0.0 0 .0 !L07 C2H3 + C2H2 C2H2 + C2H END = C4H4 + H = C4H2 + H