STR Group

Modeling of Crystal Growth and Devices

Home About Consulting Products Learn Publications Events Distributors Contact


Global Heat Transfer

     The heat transfer in an anisotropic solid body is described by the heat conductivity equation that may be written as
where T is the temperature,  is the thermal conductivity tensor and Q is the heat source due to resistive or inductive heating.
     In PVT growth, radiaton is the major mechanism of the heat transfer in gas domains. Radiative heat exchange between any solid surfaces through a non-participating fluid is modeled using the  view-factor technique on the assumption of gray-diffusive surface radiation. The total radiative flux incoming to a given surface element is calculated using the view factors . The flux incoming to the elementary surface element   is calculated as
where  is the radiation flux incident to the i-th surface element,  is the radiation flux outgoing from the i-th surface element,   is the heat emission from the i-th surface element,  is the emissivity of the i-th surface element, is the Stephan-Bolzmann constant.

Heterogeneous chemistry: generalized model

    To consider the chemical reactions running on the reactive surfaces, we have suggested an original model of the heterogeneous chemistry [1], where the total flux Ji  of the i-th species is related to its partial pressure  Pi  at a reactive surface by the Hertz-Knudsen equation while the equilibrium pressures    obey the mass action law.   is the sticking coefficient of the i-th species accounting for a possible kinetic limitation of the species adsorption/desorption rate on the reactive surface;  is the Herz-Knudsen coefficient. The model provides a unified description of the surface chemistry, applicable to various materials and growth techniques (see, for instance [2]).

        Heterogeneous processes in SiC sublimation growth

    In bulk SiC growth, the vapor includes the following species: Si, Si2C, and SiC2. The following heterogeneous reactions are assumed at the surface of SiC crystal:
An important process in SiC growth is graphitization of a SiC powder surface due to dominant Si evaporation. As a result, SiC and graphite islands coexist on the surface, in equilibrium with each other via adatom exchange. Such a surface requires consideration of the “SiC-C-vapor” three-phase equilibrium for the thermodynamic pressures     rather than of the “SiC-vapor” two-phase equilibrium corresponding to the growth of stoichiometric SiC. Then the following reaction is added to those listed above:
The sticking coefficients of all the reactive species in bulk SiC growth are assumed to be unity.

         Heterogeneous processes in AlN sublimation growth

    In bulk AlN growth, the vapor includes Al and N2.  The following heterogeneous reaction is assumed at the surface of an AlN crystal:

The surface mechanisms related to AlN growth differ significantly from those of SiC. The principal difference comes from the kinetic limitation of the adsorption/desorption rate of molecular nitrogen [3] which is a major reactive vapor species along with atomic Al. To account for this effect, we use the temperature dependent N2 sticking coefficient extracted from the data on free evaporation of AlN in vacuum [3]. This procedure provides a good agreement between the theoretical predictions and experimental data [4]. The interpretation of the specific mechanisms responsible for the surface kinetics, which is common to all Group-III nitrides, is given in [5].
    In particular, the considered kinetic factor is responsible for the extremely high growth temperatures (~1900-2300°C) conventionally used for AlN growth. The N2 sticking coefficient rises with temperature so that nitrogen becomes sufficiently reactive only at a high temperature or at a high pressure that provides a higher frequency of molecular collisions with the growth surface.

         Heterogeneous processes in GaN sublimation growth

    In bulk GaN growth, the vapor includes Ga, NH3, N2, and H2. The following heterogeneous reactions are assumed at the surface of a GaN crystal:
    The surface mechanisms related to GaN growth are similar to those of AlN. The major difference is in the number of reactive species, which requires determination of the sticking coefficients for both N2 and NH3. We use the temperature dependent N2 sticking coefficient extracted from the data on free evaporation of GaN in vacuum [5] and constant coefficient for NH3 [6].

Mass Transport in the Powder Charge

    The powder charge is considered as porous medium characterized by local porosity, granule size, and graphitization degree. Species transport in the powder is modeled in the approximation of the Darcy flow accounting for volumetric mass source due to chemical reactions on the surface of SiC or AlN granules. Temporal variation of the powder parameters (granule size, porosity and degree of SiC granule graphitization) due to granule sublimation and recrystallization is considered.
    The chemical processes on the granules are modeled using the thermodynamic model described above. In SiC growth, two types of equilibrium are considered, namely, SiC_C equilibrium accounting for the powder graphitization and SiC equilibrium corresponding to stoichiometric SiC sublimation.
The continuity equations for the whole vapor and for each species are as follows:
Flow in the porous medium is described by the Darcy law which can be written as follows:
where  is the gas viscosity, p is the gas pressure, and K is the powder permeability. The latter can be found from the granule radius and porosity using Ergun’s relationship [7]
Finally, we obtain the following equation with respect to the pressure variable:
Then the flow velocity is found from the Darcy law.

Graphitization of the SiC Powder

    The model of the powder evolution accounts for the graphitization of SiC powder granules. A graphitized granule is considered as a solid sphere consisting of the internal SiC core covered by the graphite shell (see figure below) that is permeable for the SiC vapor. Initially, the granule is not graphitized, so that the core radius  is equal to the total granule radius .
The granule surface is modeled by the three-phase SiC-C equilibrium assuming the heterogeneous reactions on both condensed phases (SiC and graphite) to run under quasi-equlibrium conditions. The solution to the mass transport problem with the SiC-C boundary conditions yields both the SiC sublimation rate  and graphite growth rate . The former quantity describes the congruent component of the granule evaporation. The latter quantity is positive if the granule evaporation is actually incongruent and the silicon atoms leave the granule at a higher rate than the carbon atoms. In the long-term powder evolution,  and   are varied independently according to the computed values of  and

    Heat exchange in the powder source and thermal insulation

One of the problems arising in modeling PVT growth is the lack of knowledge on the properties of the used materials, in particular, of the thermal insulation and powder source. The thermal insulation normally made of a graphite felt and powder are porous media with properties remarkably dependent on their microstructure [8,9]. The effective conductivity of the felt is largely determined by the heat conduction through the fibers randomly contacting with each other and through the gas. As a result, the conductivity depends much on the inert gas nature and its pressure. The electric conductivities of the powder and felt also vary with their porosity which can change during a long-term crystal growth. These factors are taken into account when modeling the temperature distribution in an inductively heated growth system.
    Another aspect, which is often neglected in simulations, is the dependence of graphite properties on its porosity. Analysis performed recently in [10] shows that the temperature dependence of the graphite heat conductivity at an arbitrary temperature can be estimated through its porosity and the conductivity measured at room temperature.

RF Source Computation

    The mathematical model for the electromagnetic field employs a quasi-stationary approach of the Maxwell equations for linear non-uniform medium. On the assumption that dependence of current intensity on time is a strictly harmonic function, the Maxwell set of equations can be written as
where is the magnetic field strength,  is the electric current density,  is electric the field intensity,    is the magnetic induction, and  is the electric induction.
    The electric current density in the overall computational domain can be expressed in terms of eddy and excitant components of the electric field as
where  is the source of the excitant electromagnetic field. The magnitude of , which is nonzero in the coil only, is determined by the inductor voltage. Finally, we can derive the following equation with respect to the field intensity

Elastic Stress Analysis

    An analysis of thermal elastic stresses in a hexagonal crystal can be carried out within the axisymmetric consideration with the symmetry axis oriented along the [0001] axis. This allows a simulation of thermal elastic stresses arising during the growth of AlN, 6H-SiC and 4H-SiC bulk crystals on the (0001) surface of the seed. The governing partial equations in cylindrical coordinates can be written as
Here, , , , and  are the stress tensor components.
    The following boundary conditions can be set for the above set of equations:
Here,  and  are the components of the unit vector normal to the surface.

    The elastic stresses are used to estimate the density of gliding dislocations in the growing crystal, which is calculated on the assumption of a full stress relaxation due to plastic deformation and crystal hardening due to dislocation formation [11]. This approach gives the following expression for the dislocation density:
where  is the applied stress,   is the critical applied stress assumed to be 106 Pa, and  is a coefficient that can be estimated by the relationship
where G is the crystal shear modulus,  is the crystal Poisson ratio, and b is the Burgers vector of a gliding dislocation.

Thermal field control by designing the crucible

    Modifying the crucible geometry is a conventional way to control the temperature distribution inside the crucible in the sublimation technique [12]. It includes searching for the optimum shape and thickness of the crucible walls, as well as artificially non-uniform thermal insulation of the crucible from the other parts of the growth system. Generally, this is a typical multi-parametric optimization problem, which is commonly solved by the trial-and-error method. Recently, an inverse-problem approach has been suggested for optimization of the crucible design [13]. This procedure has proved quite efficient, providing solutions which could hardly be found by conventional straightforward methods.
Crucible designs 1D Temperature distributions

Figure 1. Left: Conventional (left) and optimized (right) crucible designs and the respective two-dimensional temperature distributions. Right: the temperature profiles computed for the conventional (solid) and optimized (dashed) crucibles along the powder axis (top) and across the powder top (bottom)

Fig. 1 illustrates how the temperature uniformity in a SiC powder charge, necessary for its stable operation, can be improved by varying the crucible wall shape. The goal of the optimization was to provide uniform temperature distributions along the powder central axis and across its top surface. The optimal profile of the crucible wall has been obtained by the inverse-problem approach [12]. The approach can also be applied to other applications, such as (i) optimization of the thermal field during a long-term growth by moving the inductor coil and changing the inductor power supply in time, and (ii) optimization of the species mass transport at various growth stages for a long-term control of the crystal shape.


[1] A.S. Segal, A.N. Vorob’ev, S.Yu. Karpov, Yu.N. Makarov, E.N. Mokhov, M.G. Ramm, M.S. Ramm, A.D.Roenkov, Yu.A. Vodakov, and A.I. Zhmakin, Mat. Sci. and Engin. B 61-62 (1999) 40.
[2] S.Yu. Karpov, V.G. Prokofyev, E.V. Yakovlev, R.A. Talalaev, and Yu.N. Makarov, MRS J. Nitride Semicond. Res. 4 (1999) 4.
[3]L.H. Dreger, V.V. Dadape, and J.L. Margrave, J. Phys. Chem. 66 (1962) 1556.
[4] A.S. Segal, S.Yu. Karpov, Yu.N. Makarov, E.N. Mokhov, A.D. Roenkov, M.G. Ramm, Yu.A. Vodakov, J. Cryst. Growth 211 (2000) 68.
[5] M.V. Averyanova, S.Yu. Karpov, Yu.N. Makarov, I.N. Przhevalskii, M.S. Ramm, and R.A. Talalaev, MRS J. Nitride Semicond. Res. 1, (1996) 31.
[6] M. Mesrine, N. Grandjean, J. Massies, Appl. Phys. Lett. 72 (1998) 350.
[7] S. Ergun,  Chem. Engng. Progr. 48(2) (1952) 89.
[8] E.L. Kitanin, V.V. Ris, M.S. Ramm, and A.A. Schmidt, Mat. Sci. Engin. B 55 (1998) 174.
[9] E.L. Kitanin, V.V. Ris, A.A. Schmidt, S.Yu. Karpov, and M.S. Ramm, Submitted to Mat. Sci. Engin. B (2002).
[10] E.L.Kitanin, V.V. Ris, A.A. Schmidt, S.E. Demina, S.Yu. Karpov, and M.S. Ramm, ISSCRM-2002, Book of abstracts (2002) 24.
[11] . A.S. Jordan, R. Caruso, A.R.von Neida, J.W. Nielsen, J.Appl. Phys. 52 (1981) 3331.
[12] M.S. Ramm, E.N. Mokhov, S.E. Demina, M.G. Ramm, A.D. Roenkov, Yu.A. Vodakov, A.S. Segal, A.N. Vorob’ev, S.Yu. Karpov, A.V. Kulik, and Yu.N. Makarov, Mat. Sci. and Engin., B 61-62 (1999) 107.
[13] A.V. Kulik, S.E. Demina, S.K. Kochuguev, D.Kh. Ofengeim, S.Yu. Karpov, A.N. Vorob'ev, M.V. Bogdanov, M.S.Ramm, A.I. Zhmakin, A.A. Alonso, S.G. Gurevich, and Yu.N. Makarov, Mat. Res. Soc. Symp. Proc.640 (2001) H1.6.


Subscribe to STR Newsletter

Download free Software DEMO Versions and Documentation

© STR 2020. All rights reserved.