- Email: [email protected]
Transient numerical investigation of induction heating during sublimation growth of silicon carbide single crystals
Journal of Crystal Growth 247 (2003) 219–235
Transient numerical investigation of induction heating during sublimation growth of silicon carbide single crystals Olaf Klein, Peter Philip* Weierstrass Institute for Applied Analysis and Stochastics, Mohrenstrae 39, 10117 Berlin, Germany Accepted 29 August 2002 Communicated by M.E. Glicksman
Abstract This article presents transient numerical simulations of the temperature evolution during sublimation growth of SiC single crystals via physical vapor transport (also called the modified Lely method) including diffusion and radiation, investigating the influence of induction heating. Using the imposed voltage as input data, the heat sources are computed via an axisymmetric complex-valued magnetic scalar potential that is determined as the solution of an elliptic PDE. The presented results include stationary simulations of magnetic potential distributions and resulting heat sources as well as transient simulations of the temperature evolution during the heating process. We examine the effects of imposed voltage (i.e. heating power), of different coil positions, and of a moving induction coil on the evolution of the global temperature field and on the temperature at the source, at the seed, and at the blind holes. All material data used are either included or referenced. r 2002 Elsevier Science B.V. All rights reserved. PACS: 02.60.Cb; 81.10.Bk; 84.32.Hh; 44.30.+v; 44.90.+c Keywords: A1. Computer simulation; A1. Heat transfer; A1. Induction heating; growth; B2. Semiconducting silicon compounds
1. Introduction SiC bulk single crystals provide a semiconductor substrate material used in electronic and optoelectronic devices, e.g. MESFETs, thyristors, LEDs, lasers, and sensors. Due to its physical properties, SiC is especially suitable to be used in high temperature, high power, and high frequency *Corresponding author. Tel.: +49-30-20372-480; fax: +4930-204-4975. E-mail address: [email protected] (P. Philip).
A2. Growth from vapor; A2. Single crystal
applications as well as in intensive radiation environments. SiC applications require the availability of large diameter, low defect SiC boules, whose manufacturing process remains challenging inspite of progress made in recent years (cf. e.g. Ref. [1]). We consider the production of SiC single crystals by sublimation growth via physical vapor transport (PVT) (modified Lely method, see Ref. [2]). Usually, PVT growth systems consist of an induction-heated graphite crucible containing polycrystalline SiC source powder and a
0022-0248/03/$ - see front matter r 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 0 2 ) 0 1 9 0 3 - 6
220
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
single-crystalline SiC seed cooled by means of a blind hole (cf. Fig. 1). The system is degassed to 103 Pa and heated to some 1200 K to eliminate contaminants such as nitrogen. At this stage, a high-purity argon (inert gas) atmosphere is established at 105 Pa, and the temperature is further increased. At growth temperature, which can reach up to 3000 K for 6H SiC, pressure is reduced to about 2 103 Pa (cf. Ref. [3]). Due to high temperature and low pressure, SiC source powder evaporates, adding species such as Si, Si2 C; and SiC2 to the gas phase. Crystallization occurs at the cooled seed which thereby grows into the reaction chamber. As the crystal’s quality and growth rate depend strongly on the evolution of the temperature distribution (cf. Refs. [4,5]), mass transport, and species concentrations (cf. Ref. [6]), control of
these quantities is essential. However, due to the high temperatures, experimental verification of the correlation between the quantities to be regulated and control parameters such as apparatus configuration, heating power, and argon pressure is extremely difficult and costly. Thus, theoretical modeling and numerical simulation play a fundamental role in gaining understanding of the relation between control parameters and favorable growth conditions. Various recent publications report on the development and applications of numerical tools in the framework of PVT, mostly using stationary models (cf. e.g. Refs. [7–12]). Transient numerical results on the evolution of heat transfer during PVT are presented in Refs. [13,14]. Crystal growth can already occur during the heating-up stage, possibly resulting in persisting
Fig. 1. Setup of growth apparatus according to Ref. [9, Fig. 2].
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
defects such as micropipes and unwanted polytypes. Moreover, thermal stresses in the seed crystal due to temperature gradients can initiate crystal defects. Therefore one cannot restrict one’s attention to the quasi-stationary state at the end of the heating process, but it is also important to monitor and control the temperature field evolution during the heating process itself. The main objective of this work is to numerically investigate induction heating during PVT and to present systematic transient numerical parameter studies of how changes in the heating power and the coil position affect the evolution of the global temperature field and the evolution of the temperature at significant points such as the seed crystal, the SiC source, and the blind holes at the top and at the bottom of the growth apparatus (cf. Section 3.3). The numerical simulations presented in this paper are performed for the configuration depicted in Fig. 1 taken from Ref. [9, Fig. 2]. We employ the transient heat transport model previously described in detail in Refs. [14,15]. Neglecting any mechanical or chemical interactions inside both solid and gas as well as convective contributions in the gas phase, the heat transport is described by (cf. Ref. [14, (1.4), (1.5), and (3.1)]) qT divðk½b grad TÞ ¼ m½b qt ðsolid component bÞ;
r½b c½b sp
rgas
qUgas divðkgas grad TÞ ¼ 0 qt
ð1:1aÞ ðgas phaseÞ; ð1:1bÞ
where superscripts ½b refer to quantities in the solid component b (b ¼ Insulation, SiC Crystal, y), and the subscript ‘‘gas’’ refers to quantities in the gas phase. The meaning of the symbols is as follows: rgas ; r½b —mass density; c½b sp —specific heat; T—absolute temperature; t—time; kgas ; k½b —thermal conductivity; Ugas —total internal energy; m½b —power density (per volume) of heat source due to induction heating. In contrast to Ref. [14], where a homogeneous heat source was assumed inside a prescribed susceptor region, we compute the material- and temperature-dependent heat source distribution
221
according to the position of the induction coil and the voltage imposed therein (cf. Section 2). As in Ref. [14] all temperature simulations assume a gas phase consisting of pure argon. It is explained in Ref. [14, Section 5] that this is a reasonable approximation for temperature simulations, notwithstanding the fact that species such as Si, Si2 C; and SiC2 make up a significant portion of the gas mixture for high temperatures. For pure argon the internal energy is Ugas ¼
3 R T; 2 M ðArÞ
ð1:2Þ
where R denotes the universal gas constant, and M ðArÞ denotes the molecular weight of argon. Radiative heat transfer between surfaces of cavities is included using the net radiation method for diffuse-gray radiation, where the interface condition between the solid b and the gas phase reads (cf. Ref. [14, (2.1)]) ðkgas grad TÞ . ~ n gas þ Rrad Jrad ¼ ðk½b grad TÞ . ~ n gas ;
ð1:3Þ
Rrad denoting radiosity, Jrad denoting irradiation, and ~ n gas denoting the unit normal vector on the interface pointing from gas to solid. All solids are treated as opaque, except the SiC single crystal, where semi-transparency is accounted for via the band approximation model. Thus, between the SiC single crystal and adjacent solid materials, the interface condition is analogous to Eq. (1.3), whereas one has Eq. (1.3) with Rrad ¼ Jrad ¼ 0 between all other solid materials. The temperature is presumed to be continuous everywhere in the growth apparatus, disregarding the expected step temperature ðkgas grad Tgas Þ . ~ n gas ¼ xb Tgas T ½b at gas– solid interfaces, in absence of precise data for transition coefficients xb : The Stefan–Boltzmann law provides the outer boundary condition (cf. Ref. [14, (3.3)]) 4 ðkgas grad TÞ . ~ Þ; n b ¼ se½b ðTÞðT 4 Troom
ð1:4Þ
s denoting the Boltzmann radiation constant, e denoting the emissivity of the surface, and ~ nb denoting the outer unit normal vector to the solid material b: Condition (1.4) means that the growth apparatus is exposed to a black-body environment
222
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
(e.g. a large isothermal room) radiating at room temperature Troom : Analogously, we use black body phantom closures of open radiation regions such as Gtop and Gbottom depicted in Fig. 1 which emit radiation at Troom ; thereby allowing for radiative interactions between an open cavity and the ambient environment including reflections at the cavity’s surfaces. For the numerical simulations of this work, the above equations will be considered in two space dimensions, taking into account the cylindrical symmetry of the growth system.
2. The induction heating model The axisymmetric crucible is placed inside a copper induction coil as depicted in Fig. 1. An alternating voltage is imposed in the coil, resulting in an alternating current that generates a rapidly oscillating magnetic field, inducing eddy currents in the growth apparatus. The eddy currents cause heat sources due to the Joule effect. We will now describe the modeling of the induction process with the goal of computing the induced eddy currents and the resulting heat sources. It is assumed that a sinusoidal alternating voltage V ðtÞ ¼ V0 sin ðotÞ
ð2:1Þ
is imposed in the coil, since e.g. in experimental 1 setups at the Institut fur (IKZ) in . Kristallzuchtung . Berlin the voltage is more easily measured than the current. To be able to consider the problem in an axisymmetric setting, the actual coil is replaced by N cylindrical rings with a sinusoidal alternating voltage imposed in each ring. Thus, the voltage imposed in the kth ring has the form vk ðtÞ ¼ vk;0 eiot ;
ð2:2Þ
where i denotes the imaginary unit P and v1;0 ; y; vN;0 N are complex voltages satisfying k¼1 vk;0 ¼ V0 : Here, as in the following, we use the complex representation of sinusoidal functions. All solid materials in the growth system are considered as potential conductors, whereas the 1
Institute of Crystal Growth.
gas phase is treated as a perfect insulator. We assume that displacement currents as well as surface currents can be neglected. All the involved electromagnetic quantities are assumed to be axisymmetric and sinusoidal. Using cylindrical coordinates ðr; W; zÞ and their standard basis e W ;~ e z Þ; the complex current density is supposed ð~ e r ;~ to be given by ~ j ¼ jðr; zÞ~ e W : It is shown in Refs. [16,17] that under these assumptions, Maxwell’s equations and Ohm’s law allow one to find a complex-valued magnetic scalar potential fðr; zÞ ~ ¼ eiot fðr; zÞ~ ~ is the magsuch that A e W ; where A netic vector potential, i.e. the magnetic induction ~: ~ ¼ curl A reads B The potential f is determined from the system of elliptic partial differential equations [17, (22), (29)], which we rewrite in the following divergence form (2.3) which is more suitable for our numerical approach via a finite volume discretization: gradðrfÞ n div ¼ 0 in the gas phase; ð2:3aÞ r2 gradðrfÞ i osc f sc vk;0 ¼ þ r2 r 2pr2 in the kth coil ring;
ð2:3bÞ
gradðrfÞ i osc f ¼0 þ r2 r in other conducting materials;
ð2:3cÞ
n div
n div
where n denotes the reciprocal of the magnetic permeability, and sc is the electrical conductivity. The quantities n and sc can vary in space, but they are supposed to be constant in time. System (2.3) is completed by interface and boundary conditions. Owing to the assumption of no surface currents, we have the interface conditions [Ref. 17, Eq. (30)]: npMaterial1 .~ gradðrfÞp n Material1 Material1 r2 npMaterial2 .~ ¼ gradðrfÞp n Material1 Material2 r2 ð2:4Þ on interfaces between Material1 and Material2 ; where p denotes the restriction to the respective material, and ~ n Material1 denotes the outer unit
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
normal vector to Material1 : It is also assumed that f is continuous throughout the whole domain and that f ¼ 0 both on the symmetry axis r ¼ 0 and sufficiently far from the growth apparatus (in Section 3.2, it is found that using a domain radius and height of 1.2 and 1.8 m, respectively, is sufficiently accurate for our purposes). From Ref. [17, Eq. (28)], we get ~ j ¼ j0 eiot~ eW; with j0 ¼ 8 sc vk;0 > < io sc f þ 2pr io sc f > :
in the kth coil ring; in other conducting
ð2:5Þ
materials:
Let Ok denote the two-dimensional domain of the kth coil ring. The total current in the kth ring is then given by eiot Jk ; with Z j0 ðxÞ dx: ð2:6Þ Jk :¼ Ok
Since the coil rings constitute the approximation of a single, connected coil, the total current must be the same in each coil ring, i.e. Jk ¼ Jkþ1
for
k ¼ 1; y; N 1:
ð2:7Þ
Hence, we have to find a complex decomposition V0 ; such that v1;0 ; y; vN;0 of the imposed voltage P Eq. (2.7) is satisfied as well as N v k¼1 k;0 ¼ V0 : This leads to a linear system for v1;0 ; y; vk;0 (see Ref. [18] for more details). After solving the combined system for ðf; v1;0 ; y; vN;0 Þ; we use Eq. (2.5) to compute the power density of the heat sources according to mðt; r; zÞ ¼
jj0 ðr; zÞj2 : 2 sc ðTðt; r; zÞÞ
3. Numerical experiments 3.1. General setup and methods All numerical simulations presented in the following were performed for a growth system consisting of a container with a height of 25 cm and a radius of 8.4 cm placed inside of 5 hollow rectangular-shaped copper induction rings as displayed in Fig. 1, each ring having a surface thickness of 1 mm. The geometric proportions of the coil rings are depicted in Fig. 2. While the horizontal gap between container and coil is fixed at 4.2 cm in all our experiments, the vertical coil position is varied for the simulations in Section 3.3. The prescribed heating voltage satisfies Eq. (2.1). We fix its frequency at f ¼ 10 kHz, whereas we vary the effective voltage such that 190 VpVeff p270 V: The angular frequency and the p amplitude are provided by o ¼ 2pf and V0 ¼ ffiffiffi Veff 2: It is reiterated that only argon is considered in the gas phase and that it is the only material considered a perfect insulator. The solid materials are assumed to be homogeneous and pure materials. In particular, inhomogeneities in the source powder, and chemical changes occurring during the growth process such as sintering and
ð2:8Þ
The assumption of a sinusoidal form of the electromagnetic quantities is only admissible if the involved material parameters are constant in time, and if there is no movement. Thus, to account for the temperature-dependence of the electrical conductivity and for changing coil positions during the transient simulations in Section 3.3, the quasi-stationary electromagnetic problem presented above, including the decomposition of the imposed voltages, is solved in each time step of the transient computations for the heat evolution problem.
223
Fig. 2. Geometric proportions of induction coil rings.
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
224
graphitization of the source, and Si accumulating in the insulation felt are ignored. In real growth systems the coil is usually cooled very efficiently, e.g. by water flowing inside the coil rings. Hence it is reasonable to assume that the coil is always at room temperature, its material data not being changed during the experiments. The complete list of material data used for our numerical simulations of induction heating can be found in Appendix A. The material parameters used for the simulations of the temperature field evolution are identical to those used in Ref. [14] and can be found in Ref. [14, Appendix A]. A finite volume method based on Ref. [19] is used for the space discretizations arising in the stationary computation of the magnetic scalar potential and in the transient temperature simulations, where an implicit Euler scheme provides the time discretization. The discrete scheme has been implemented using the pdelib program package, being developed at the Weierstrass Institute of Applied Analysis and Stochastics (WIAS), Berlin (cf. Ref. [20]). For the numerical methods used for radiation and view factor calculation we refer to Ref. [14, Section 4]. 3.2. Stationary computations of heat sources This section presents and discusses results of three stationary computations at Troom ¼ 293 K of the magnetic scalar potential f and the resulting power densities m: In addition, the total heating ½apparatus powers Ptotal and P½coil total are computed in the growth apparatus and in the coil rings, respectively, by integrating m over the corresponding regions. More precisely Z X ½apparatus ¼ m½b ; Ptotal P½coil total ¼
Z
solids b a coil rings
m½Cu :
solid b
ð3:1Þ
coil rings
In practice the integrals in Eq. (3.1) are approximated by sums over finite volume elements. The simulations use Veff ¼ 230 V, the coil’s upper rim is at z ¼ 14 cm and its lower rim is level with the lower rim of the growth apparatus at z ¼ 0: The results show that even though the
outcome for the magnetic potential distribution does depend on the size of the domain used in the computation, the arising heat source distribution is not affected significantly. The simulations were performed on domains Osmall (radius rsmall ¼ 0:4 m; height hsmall ¼ 1:0 m), Omedium (rmedium ¼ 1:2 m; hmedium ¼ 1:8 m), and Olarge (rlarge ¼ 2:8 m; hlarge ¼ 3:4 m), the growth apparatus being always approximately in the center of the domain. It is noted that instead of imposing the boundary condition f ¼ 0 at r ¼ 0; f ¼ 0 is set at the vertical line r ¼ 1 mm, where f is considered to be sufficiently small. This permits the evaluation of terms with r in the denominator (cf. Eqs. (2.3) and (2.4)). Fig. 3 depicts the results of the three numerical experiments, Osmall in column (a), Omedium in column (b), and Olarge in column (c). Row (1) shows the imaginary part ImðfÞ of the magnetic scalar potential with isolevels spaced at Diso ¼ 105 Wb/m, row (2) shows the real part ReðfÞ (Diso ¼ 107 Wb/m), and row (3) shows the induced power densities m (Diso ¼ 2 MW=m3 ). The pictures demonstrate that ImðfÞ reaches rather far into the surrounding space of the growth apparatus, and the calculated global distribution depends noticeably on the size of the domain. The reach of ReðfÞ is considerably less, as can be seen from the agreement of (2)(b) and (2)(c) in Fig. 3. However, in a close neighborhood of the growth apparatus and the coil rings the distributions of ImðfÞ and ReðfÞ are almost identical for the three domains. This results in virtually identical power density distributions (see row (3) of Fig. 3). Moreover, the values in Table 1 show that the ranges of ImðfÞ; ReðfÞ; and m (both in the coil rings and in the apparatus) agree up to 3 digits between Omedium and Olarge : The same holds true ½apparatus for the total heating powers P½coil : total and Ptotal In consequence Omedium provides a reasonable compromise between accuracy and computing expenditure and is thus used for all computations of heat sources in Section 3.3. The results depicted in Fig. 3 concur qualitatively with the results found in Refs. [9, Fig. 3, 13, Fig. 2]. While a frequency of 10 kHz is not high enough to produce a full-fledged skin effect, the heat
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
225
Fig. 3. Magnetic scalar potential f (1st row: imaginary part, 2nd row: real part) and resulting power density m of heat sources (3rd row) for three different domain sizes. See text for discussion of results.
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
226
Table 1 Extreme values of the real and imaginary parts of the magnetic scalar potential, maximal values of power densities, and values of total powers established during numerical simulations on Osmall ; Omedium ; and Olarge
6
ReðfÞmin ð10 Wb=mÞ ReðfÞmax ð106 Wb=mÞ ImðfÞmin ð104 Wb=mÞ ImðfÞmax ð104 Wb=mÞ 8 3 m½coil max ð10 W=m Þ 7 3 ð10 W=m Þ m½apparatus max P½coil ðkWÞ total P½apparatus ðkWÞ total ½apparatus P½coil þ P ðkWÞ total total
Osmall
Omedium
Olarge
8.80 4.73 1.79 0.00 3.49 1.45 1.63 4.83 6.46
8.82 4.68 1.78 0.00 3.39 1.46 1.57 4.89 6.46
8.82 4.68 1.78 0.00 3.39 1.46 1.57 4.89 6.46
sources decay quickly below a conductor’s surface (see Fig. 3(3), where the maximal power density is established in the lower right-hand corner of the graphite). ½apparatus The values for P½coil in Table 1 total and Ptotal ½apparatus show that some 25 percent of P½coil is total þ Ptotal lost in the coil rings. 3.3. Transient simulation of the temperature evolution In this section we report on transient numerical simulations of the temperature evolution in the growth apparatus, each run starting at T ¼ 293 K. We varied the parameters for argon pressure, heating voltage, and vertical coil positions, to assess their influence during the heating process. It was found in Section 3.2 that the domain Omedium of radius 1.2 m and height 1.8 m is suitable for accurate heat source computations. Hence Omedium is used for all heat source simulations in the current section. As already pointed out at the end of Section 2, the heat source distribution is computed in each time step of the transient simulations to allow for the temperature dependence of electrical conductivities (cf. Eqs. (A.1)). The temperature distribution is only calculated in the growth apparatus, i.e. on a much smaller grid GT than the grid Gm used for the heat source computations. In consequence in each time step data for T are carried from GT to Gm and data for m are carried from Gm to GT :
It is noted that the relaxation times of the pressure and of the temperature in the gas phase always lie orders of magnitude below the relaxation times of the temperature in the solid components. Since the gas does not have a significant influence on the temperature of the solid parts, each temporal snapshot of a transient simulation furnishes the quasi-stationary temperature distribution in the gas phase determined by the temperature distribution on the adjacent solid surfaces. In particular, for a constant pressure pðArÞ ; the time derivative in Eq. (1.1b) vanishes identically as can be seen from Eq. (1.2) and the gas law pðArÞ ¼ rðArÞ ðR=M ðArÞ ÞT: Thus, the quasistationary temperature distribution in the gas phase does not depend on the gas pressure. The argument of the preceding paragraph suggests that Eq. (1.1a) could be coupled with the stationary version of Eq. (1.1b), i.e. with divðkðArÞ grad TÞ ¼ 0: However, for the numerical regularity of the problem, it is advantageous to keep the time derivative in the gas equation which is also helpful within the present code implementation. For the numerical experiments described in the following, instead of using a given argon pressure, we fixed rðArÞ ¼ 3:73 103 kg=m3 ; which corresponds to T ¼ 2575 K and pðArÞ ¼ 2 103 Pa. The total heating power in the growth apparatus P½apparatus and in the coil rings P½coil total total changes during each simulation run due to the temperature-dependent electrical conductivities. Moreover, if present, a movement of the induction coil also changes P½apparatus and P½coil total total : For one of our simulations, using Veff ¼ 230 V and the fixed coil position C16 2 according to the notation introduced in the next paragraph, Figs. 4(a)(1) and (a)(2) depict the time evolution of P½apparatus and P½coil total total , respectively. The curves show an inverse behavior, i.e. P½apparatus decreases as P½coil total total increases and vice versa. The peaks after some 500 s in Fig. 4(a) correspond to the peaks in Fig. 4(b) which show the electrical conductivity of the graphite crucible depending on temperature in Fig. 4(b)(1) and depending on the maximal temperature established in the graphite at time t in Fig. 4(b)(2). The minimum in Fig. 4(a)(1) matches the maximum in Fig. 4(b)(2), showing that the system is in a state
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
227
Fig. 4. Time evolution of total heating power in growth apparatus (a)(1) and in coil rings (a)(2) for numerical experiment C16 2 230 V in comparison with the temperature-dependent electrical conductivity of the graphite crucible (b).
where increasing the electrical conductivity of the apparatus (lowering its resistance) reduces ½apparatus Ptotal : We now discuss three series of numerical 16 18 experiments, referred to as C14 0 ; C2 ; and C4 ; respectively, each series using a different vertical coil position. For C14 the coil is positioned 0 between z ¼ 0 and 14 cm, i.e. the lower rim of the bottom coil ring is at z ¼ 0; and the upper rim of the top coil ring is at z ¼ 14 cm. The meaning of 18 C16 2 and C4 is analogous. There are five experiments in each series using different values for Veff ; namely Veff ¼ 190, 210, 230, 250, and 270 V. Subsequently, a fourth series C18-14 4-0 is considered, using a moving induction coil. For each simulation, Table 2 contains P½apparatus and P½coil total total ; i.e. the
total heating power in the growth apparatus and in the coil rings, respectively, established in the stationary final state. Results for the time evolution and the final values of the temperature evolution at four points of particular importance are depicted in Figs. 5–7, respectively. The monitored points are located at the blind holes at the top and at the bottom of the growth apparatus, and at the surface of the crystal seed and of the source powder, respectively (labeled Ttop ; Tbottom ; Tseed ; and Tsource in Fig. 1). The significance of Ttop and Tbottom lies in the blind holes being the principal locations for temperature measurements; Tseed and Tsource are of interest, since their difference is a key control factor for the crystal’s growth rate and quality (cf. e.g.
228
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
Fig. 5. Time evolution of Ttop (row (1)), Tbottom (row (2)), Tseed (row (3)), and Tsource Tseed (row (4)) for the simulation series C14 0 (column (a)) and C16 2 (column (b)).
Refs. [2,5]). Results of the time evolution of the global temperature field and the heat sources are shown in Fig. 8, and results of the quasi-stationary
temperature distributions at the crystal and its immediate surroundings are depicted in Fig. 9. We now proceed to discuss our results in detail.
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
229
Fig. 6. Time evolution of Ttop (row (1)), Tbottom (row (2)), Tseed (row (3)), and Tsource Tseed (row (4)) for simulation series C18 4 (column 18-14 (a)) and C4-10 (column (b)).
We find that for each time t; the temperatures Ttop ; Tbottom ; Tseed ; and Tsource depend nearly linearly on the heating voltage Veff : For Ttop ;
Tbottom ; and Tseed ; this can be observed in Figs. 5 and 6. Moreover, for the quasi-stationary state at t ¼ 30 000 s, the nearly linear dependence of these
230
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
Fig. 7. Temperatures Ttop (in (1)), Tbottom (in (2)), Tseed (in (3)), and Tsource Tseed (in (4)) for the quasi-stationary state at t ¼ 30 000 s 16 18 for the simulation series C14 0 ; C2 ; and C4 depending on the effective heating voltage.
temperatures and also of Tsource Tseed is illustrated in Fig. 7. Figs. 5–7 show furthermore that Ttop ; Tseed ; and Tsource increase significantly be18 tween C14 0 (low coil position) and C4 (high coil position), whereas Tbottom remains almost constant. This can be explained by the heat sources being shifted upwards and the source powder’s insulating property which, though less prominent, is still effective even at high temperatures (cf. Fig. 8(a) and Ref. [14, (A.6b)]), thereby hindering heat generated below the source powder from reaching regions above the source powder.
The shape of the curves in row (4) of Figs. 5 and 6 depicting Tsource Tseed is also due to the insulating characteristic of the source powder and the gas phase, initially keeping the powder’s surface below the temperature of the seed which the heat reaches via conduction through the graphite. This causes the negative peak in the corresponding graphs. At higher temperatures, radiative heat transfer becomes more effective both in the gas phase and the porous source, resulting in the source growing warmer than the seed, as is required for crystal growth. Moreover,
18-14 Fig. 8. Time evolution of heating process for numerical experiment C4-0 230 V (coil moving at 1:33 cm/h). Column (a): temperature evolution, difference between neighboring isotherms is 20 K. Column (b): heat source evolution, difference between isolevels is 10 kW/m3 ; darker regions indicate larger heat sources.
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
231
232
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
Fig. 9. Close-ups of the temperature distribution in the crystal, the gas phase, and the crucible directly above the crystal for the quasistationary states at t ¼ 30 000 s of the numerical experiments with lowest and medium coil position in (a) and (b), respectively. Temperature difference between neighboring isotherms is 5 K.
it is seen in Fig. 7(4) that at the quasi-stationary state, a lower coil position or a higher heating voltage result in Tsource Tseed being increased as heat reaches the powder’s surface more readily at higher temperatures due to radiative heat transfer. Comparing our findings in Fig. 7(4) with [21, Fig. 6] establishes a qualitative agreement, whereas a quantitative comparison is meaningless due to the different setups. Fig. 7(4) indicates that Tsource Tseed can be tuned effectively by controling the heating voltage and the vertical coil position. We now discuss the series of numerical experiments C18-14 4-0 ; where the coil starts at the position used for C18 4 ; moves downwards at the rate 1.33 cm/h, and stops at the position used for C14 0 : During these simulations the grid for the heat source computations is newly generated in each time step. The objective of the series C18-14 is to 4-0 study how a moving coil affects the temperature field evolution during the heating process. In this article, we do not investigate coil movements that follow the growing crystal in the quasi-stationary temperature field which are often used in growth experiments to compensate the influence of the growing crystal on the temperature distribution. The rate of such coil movements is usually in the order of 1 mm/h. As expected, the curves for the moving coil (Fig. 6(b)) initially coincide with the corresponding curves of C18 4 (Fig. 6(a)) and coincide with the corresponding curves of C14 0 (Fig. 5(a)) after the coil’s lowest position is reached at t ¼ 3 h: For experiment C18-14 4-0 230 V, we discuss our results concerning the time evolution of the global temperature distribution and the heat sources. It can be seen in column (a) of Fig. 8 that the minimal temperature Tmin is always established at
the outside of the outer insulation material, where the isotherms become very dense at higher temperatures, producing the dark outer regions in Fig. 8(2)(a) and (3)(a). The maximal temperature Tmax is found in the graphite strip between source and insulation in Fig. 8(1)(a) and inside the labeled isotherms below the source in Fig. 8(2)(a) and (3)(a). The change of the location of Tmax corresponds to the heat sources moving downwards together with the induction coil according to Fig. 8(b). Fig. 8(a) also demonstrates how the insulating property of the gas phase and the source powder changes with time, causing a large temperature gradient in the gas phase in Fig. 8(1)(a) and local minima in the source in Figs. 8(1)(a) and (2)(a). Due to radiative heat transfer through the gas phase and through the pores of the source, the temperature on the boundary of the growth chamber (and hence also inside the gas itself) is more homogeneous in Figs. 8(2)(a) and (3)(a). Moreover, the local minimum has vanished in Fig. 8(3)(a), even though some of the source’s insulating quality persists. The increase in P½apparatus between Figs. 8(1)(b) total and (2)(b) followed by the drop between Figs. 8(2)(b) and (3)(b) is consistent with Fig. 4(a)(1) and Table 2: The increase corresponds to the rise after the initial minimum in Fig. 4(a)(1), and the drop matches the power’s decrease with a lower coil position found in Table 2. Finally, Fig. 9 displays close-ups of the final temperature distributions in the crystal’s immediate surroundings for two coil positions 16 C14 0 230 V (in Fig. 9(a)) and C2 230 V (in Fig. 9(b)). The picture in Fig. 9(a) also constitutes a close-up of Fig. 8(3)(a), since the stationary states are identical for the simulations C14 0 230 V
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
233
Table 2 Total heating power P½apparatus in the growth apparatus and P½coil total total in the coil rings for the quasi-stationary state at t ¼ 30 000 s, depending on heating voltage and coil positions Series
Veff P½apparatus ðkWÞ total
C14 0 C16 2 C18 4 C18-14 4-0
P½coil total ðkWÞ
190
210
230
250
270
190
210
230
250
270
3.15 3.40 3.55 3.15
3.87 4.19 4.37 3.87
4.67 5.06 5.27 4.67
5.55 6.01 6.26 5.55
6.51 7.05 7.34 6.51
1.08 1.09 1.10 1.08
1.32 1.33 1.37 1.32
1.58 1.59 1.60 1.58
1.86 1.88 1.89 1.86
2.17 2.18 2.20 2.17
and C18-14 4-0 230 V. The lower temperatures found for the lower coil position in Fig. 9, are once more explained by the powder source forming a barrier for heat generated in the lower part of the apparatus. Moreover, Fig. 9 shows that the final temperature gradient between source and seed is nearly independent of the coil position and approximately constant along the r ¼ 0 axis with a slight increase at the seed.
. Bottcher, Thomas Muller, . Detlev Schulz, and Dietmar Siche. We also thank Jurgen . Sprekels and Manfred Uhle of the WIAS in Berlin for helpful discussions and advice. This research was partly funded by the Bundesministerium fur . Bildung, Wissenschaft, Forschung und Technologie2 (BMBF) within the programs Neue Mathematische Verfahren in Industrie und Dienstleistungen3 # 03SPM3B5 and Mathematische Verfahren zur Losung von Problemstellungen in . Industrie und Wirtschaft 4 # 03SP7FV16.
4. Conclusions The influence of induction heating during PVT was studied by means of transient numerical simulations including diffusive and radiative heat transfer, where a quasi-stationary elliptic problem is solved in each time step to determine temperature- and coil-position-dependent heat sources. We presented results of computed temperature evolutions depending on heating voltage and vertical coil position. The results agree with physical expectations and provide insights on how to control both the global temperature distribution evolution and the temperature evolution at the seed and the source by adjusting heating voltage and vertical coil position.
Appendix A. Material data for induction heating
Acknowledgements
2 German Ministry for Education, Science, Research, and Technology. 3 New Mathematical Methods in Manufacturing- and Service Industry. 4 Mathematical Methods Solving Problems in Industry and Business.
The IKZ in Berlin has provided us with material data and experimental data throughout the project. In particular, we wish to thank Klaus
The relevant materials are argon, which is the only material treated as a perfect insulator, copper, graphite crucible, insulation, SiC powder, and SiC crystal. In each material the reciprocal of the magnetic permeability is assumed to be n ¼ m1 0 ¼ ð107 =4pÞ A m=V s: A.1. Mass densities of conducting materials The values of Table 3 are according to Refs. [22 (Cu); 23 (graphite crucible, insulation); 24 (SiC powder); 25 (SiC crystal)].
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
234
Table 3 Mass densities of conducting materials Material
Cu
Graphite crucible
Insulation
SiC powder
SiC crystal
r ½kg=m3
8930
1750
170
1700
3140
A.2. Electrical conductivities of conducting materials According to Ref. [22] the value for the electrical conductivity of copper at room temperature is s½Cu ¼ ð1:7 108 O mÞ1 ¼ 5:9 107 ð1=O mÞ: The material parameters used for the electrical conductivities of the remaining conducting materials are s½Crucible ðTÞ ¼ !!1 6 ln ðT=1023 KÞ 2 10 ; 28:9 18:8 exp 2:37 Om ðA:1aÞ s½Insulation ðTÞ ¼
245:7 1 ; 1 þ T=2500 K O m
s½SiCPowder ðTÞ ¼ 100
1 ; Om
ðA:1bÞ ðA:1cÞ
1 : ðA:1dÞ Om Eq. (A.1a) is according to Ref. [26], and Eqs. (A.1b), (A.1c), and (A.1d) are according to Ref. [27]. s½SiCCrystal ðTÞ ¼ 105
References [1] St.G. Muller, . R.C. Glass, H.M. Hobgood, V.F. Tsvetkov, M. Brady, D. Henshall, J.R. Jenny, D. Malta, C.H. Carter Jr., J. Crystal Growth 211 (2000) 325. [2] A.O. Konstantinov, Sublimation growth of SiC, in: G.L. Harris (Ed.), Properties of Silicon Carbide, EMIS Datareview Series, no. 13, Institution of Electrical Engineers, INSPEC, London, pp. 170–203. [3] D.L. Barrett, J.P. McHugh, H.M. Hobgood, R.H. Hopkins, P.G. McMullin, R.C. Clarke, W.J. Choyke, J. Crystal Growth 128 (1993) 358. . [4] H.-J. Rost, D. Siche, J. Dolle, W. Eiserbeck, T. Muller, D. Schulz, G. Wagner, J. Wollweber, Mater. Sci. Eng. B 61–62 (1999) 68.
[5] N. Schulze, D.L. Barrett, G. Pensl, Appl. Phys. Lett. 72 (13) (1998) 1632. [6] A.S. Segal, A.N. Vorob’ev, S.Yu. Karpov, E.N. Mokhov, M.G. Ramm, M.S. Ramm, A.D. Roenkov, Yu.A. Vodakov, Yu.N. Makarov, J. Crystal Growth 208 (2000) 431. [7] S.Yu. Karpov, Yu.N. Makarov, M.G. Ramm, Phys. Stat. Sol. B 202 (1997) 201. [8] K. Chourou, M. Anikin, J.M. Bluet, J.M. Dedulle, R. Madar, M. Pons, E. Blanquet, C. Bernard, P. Grosse, C. Faure, G. Basset, Y. Grange, Mater. Sci. Eng. B 61–62 (1999) 82. [9] M. Pons, M. Anikin, K. Chourou, J.M. Dedulle, R. Madar, E. Blanquet, A. Pisch, C. Bernard, P. Grosse, C. Faure, G. Basset, Y. Grange, Mater. Sci. Eng. B 61–62 (1999) 18. [10] 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, Yu.N. Makarov, Mater. Sci. Eng. B 61–62 (1999) 107. [11] S.Yu. Karpov, A.V. Kulik, I.A. Zhmakin, Yu.N. Makarov, E.N. Mokhov, M.G. Ramm, M.S. Ramm, A.D. Roenkov, Yu.A. Vodakov, J. Crystal Growth 211 (2000) 347. [12] M. Selder, L. Kadinski, Yu. Makarov, F. Durst, P. Wellmann, T. Straubinger, D. Hoffmann, S. Karpov, M. Ramm, J. Crystal Growth 211 (2000) 333. [13] Q.-S. Chen, H. Zhang, V. Prasad, C.M. Balkas, N.K. Yushin, A System Model for Silicon Carbide Crystal Growth by Physical Vapor Transport Method, 1999 National Heat Transfer Conference, NHTC99-222, ASME, 1999, pp. 1–8. ! [14] O. Klein, P. Philip, J. Sprekels, K. Wilmanski, J. Crystal Growth 222 (2001) 832. ! [15] N. Bubner, O. Klein, P. Philip, J. Sprekels, K. Wilmanski, J. Crystal Growth 205 (1999) 294. [16] S. Clain, J. Rappaz, M. Swierkosz, Coupling between nonlinear Maxwell and heat equations for an induction heating problem: modeling and numerical methods, in: M. Krizek et al. (Ed.), Finite element methods. 50 years of the Courant element, Lecture Notes in Pure and Applied Mathematics Vol. 164, Marcel Dekker, New York, NY, pp. 163–171. [17] J. Rappaz, M. Swierkosz, Modelling in numerical simulation of electromagnetic heating, in: Modelling and Optimization of Distributed Parameter Systems, Warsaw, 1995, Chapman & Hall, New York, pp. 313–320. [18] O. Klein, P. Philip, IEEE Trans. Mag. 38 (3) (2002) 1519. [19] J. Fuhrmann, On numerical solution methods for nonlinear parabolic problems, in: Modeling and Computation
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235 in Environmental Sciences, Notes on Numerical Fluid Mechanics, Vol. 59, Vieweg, Braunschweig, Germany, pp. 170–180. [20] J. Fuhrmann, Th. Koprucki, H. Langmach, pdelib: an open modular tool box for the numerical solution of partial differential equations. Design patterns, Proceedings of the 14th GAMM Seminar on Concepts of Numerical Software, Kiel, January 23–25, 1998 (Kiel, Germany), 2001. [21] D.L. Barrett, R.G. Seidensticker, W. Gaida, R.H. Hopkins, W.J. Choyke, J. Crystal Growth 109 (1991) 17. [22] G. Benkowsky, Induktionserw.armung. H.arten, Gluhen, . . Schmelzen, Loten, Schweien, 5th edition, Verlag Technik GmbH, Berlin, 1990.
235
[23] T. Muller, . D. Schulz, D. Siche, Material vendor’s data, communicated by the Institut fur . Kristallzuchtung, . Berlin, 1999. [24] T. Muller, . D. Schulz, D. Siche, Institut fur . Kristallzuchtung, . Berlin, 1999, unpublished experimental data. [25] O. Nilsson, H. Mehling, R. Horn, J. Fricke, R. Hofmann, S.G. Muller, . R. Eckstein, D. Hofmann, High Temp. – High Press. 29 (1) (1997) 73. [26] J.G. Hust (Ed.), Standard Reference materials: finegrained, isotropic graphite for use as NBS thermophysical property RM’s from 5 to 2500 K, Tech. Report NBS/SP260/89, National Bureau of Standards, September 1984. [27] P. R(aback, Modeling of the sublimation growth of silicon carbide crystals, Ph.D. Thesis, Helsinki University of Technology, 1996.
Transient numerical investigation of induction heating during sublimation growth of silicon carbide single crystals Olaf Klein, Peter Philip* Weierstrass Institute for Applied Analysis and Stochastics, Mohrenstrae 39, 10117 Berlin, Germany Accepted 29 August 2002 Communicated by M.E. Glicksman
Abstract This article presents transient numerical simulations of the temperature evolution during sublimation growth of SiC single crystals via physical vapor transport (also called the modified Lely method) including diffusion and radiation, investigating the influence of induction heating. Using the imposed voltage as input data, the heat sources are computed via an axisymmetric complex-valued magnetic scalar potential that is determined as the solution of an elliptic PDE. The presented results include stationary simulations of magnetic potential distributions and resulting heat sources as well as transient simulations of the temperature evolution during the heating process. We examine the effects of imposed voltage (i.e. heating power), of different coil positions, and of a moving induction coil on the evolution of the global temperature field and on the temperature at the source, at the seed, and at the blind holes. All material data used are either included or referenced. r 2002 Elsevier Science B.V. All rights reserved. PACS: 02.60.Cb; 81.10.Bk; 84.32.Hh; 44.30.+v; 44.90.+c Keywords: A1. Computer simulation; A1. Heat transfer; A1. Induction heating; growth; B2. Semiconducting silicon compounds
1. Introduction SiC bulk single crystals provide a semiconductor substrate material used in electronic and optoelectronic devices, e.g. MESFETs, thyristors, LEDs, lasers, and sensors. Due to its physical properties, SiC is especially suitable to be used in high temperature, high power, and high frequency *Corresponding author. Tel.: +49-30-20372-480; fax: +4930-204-4975. E-mail address: [email protected] (P. Philip).
A2. Growth from vapor; A2. Single crystal
applications as well as in intensive radiation environments. SiC applications require the availability of large diameter, low defect SiC boules, whose manufacturing process remains challenging inspite of progress made in recent years (cf. e.g. Ref. [1]). We consider the production of SiC single crystals by sublimation growth via physical vapor transport (PVT) (modified Lely method, see Ref. [2]). Usually, PVT growth systems consist of an induction-heated graphite crucible containing polycrystalline SiC source powder and a
0022-0248/03/$ - see front matter r 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 0 2 ) 0 1 9 0 3 - 6
220
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
single-crystalline SiC seed cooled by means of a blind hole (cf. Fig. 1). The system is degassed to 103 Pa and heated to some 1200 K to eliminate contaminants such as nitrogen. At this stage, a high-purity argon (inert gas) atmosphere is established at 105 Pa, and the temperature is further increased. At growth temperature, which can reach up to 3000 K for 6H SiC, pressure is reduced to about 2 103 Pa (cf. Ref. [3]). Due to high temperature and low pressure, SiC source powder evaporates, adding species such as Si, Si2 C; and SiC2 to the gas phase. Crystallization occurs at the cooled seed which thereby grows into the reaction chamber. As the crystal’s quality and growth rate depend strongly on the evolution of the temperature distribution (cf. Refs. [4,5]), mass transport, and species concentrations (cf. Ref. [6]), control of
these quantities is essential. However, due to the high temperatures, experimental verification of the correlation between the quantities to be regulated and control parameters such as apparatus configuration, heating power, and argon pressure is extremely difficult and costly. Thus, theoretical modeling and numerical simulation play a fundamental role in gaining understanding of the relation between control parameters and favorable growth conditions. Various recent publications report on the development and applications of numerical tools in the framework of PVT, mostly using stationary models (cf. e.g. Refs. [7–12]). Transient numerical results on the evolution of heat transfer during PVT are presented in Refs. [13,14]. Crystal growth can already occur during the heating-up stage, possibly resulting in persisting
Fig. 1. Setup of growth apparatus according to Ref. [9, Fig. 2].
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
defects such as micropipes and unwanted polytypes. Moreover, thermal stresses in the seed crystal due to temperature gradients can initiate crystal defects. Therefore one cannot restrict one’s attention to the quasi-stationary state at the end of the heating process, but it is also important to monitor and control the temperature field evolution during the heating process itself. The main objective of this work is to numerically investigate induction heating during PVT and to present systematic transient numerical parameter studies of how changes in the heating power and the coil position affect the evolution of the global temperature field and the evolution of the temperature at significant points such as the seed crystal, the SiC source, and the blind holes at the top and at the bottom of the growth apparatus (cf. Section 3.3). The numerical simulations presented in this paper are performed for the configuration depicted in Fig. 1 taken from Ref. [9, Fig. 2]. We employ the transient heat transport model previously described in detail in Refs. [14,15]. Neglecting any mechanical or chemical interactions inside both solid and gas as well as convective contributions in the gas phase, the heat transport is described by (cf. Ref. [14, (1.4), (1.5), and (3.1)]) qT divðk½b grad TÞ ¼ m½b qt ðsolid component bÞ;
r½b c½b sp
rgas
qUgas divðkgas grad TÞ ¼ 0 qt
ð1:1aÞ ðgas phaseÞ; ð1:1bÞ
where superscripts ½b refer to quantities in the solid component b (b ¼ Insulation, SiC Crystal, y), and the subscript ‘‘gas’’ refers to quantities in the gas phase. The meaning of the symbols is as follows: rgas ; r½b —mass density; c½b sp —specific heat; T—absolute temperature; t—time; kgas ; k½b —thermal conductivity; Ugas —total internal energy; m½b —power density (per volume) of heat source due to induction heating. In contrast to Ref. [14], where a homogeneous heat source was assumed inside a prescribed susceptor region, we compute the material- and temperature-dependent heat source distribution
221
according to the position of the induction coil and the voltage imposed therein (cf. Section 2). As in Ref. [14] all temperature simulations assume a gas phase consisting of pure argon. It is explained in Ref. [14, Section 5] that this is a reasonable approximation for temperature simulations, notwithstanding the fact that species such as Si, Si2 C; and SiC2 make up a significant portion of the gas mixture for high temperatures. For pure argon the internal energy is Ugas ¼
3 R T; 2 M ðArÞ
ð1:2Þ
where R denotes the universal gas constant, and M ðArÞ denotes the molecular weight of argon. Radiative heat transfer between surfaces of cavities is included using the net radiation method for diffuse-gray radiation, where the interface condition between the solid b and the gas phase reads (cf. Ref. [14, (2.1)]) ðkgas grad TÞ . ~ n gas þ Rrad Jrad ¼ ðk½b grad TÞ . ~ n gas ;
ð1:3Þ
Rrad denoting radiosity, Jrad denoting irradiation, and ~ n gas denoting the unit normal vector on the interface pointing from gas to solid. All solids are treated as opaque, except the SiC single crystal, where semi-transparency is accounted for via the band approximation model. Thus, between the SiC single crystal and adjacent solid materials, the interface condition is analogous to Eq. (1.3), whereas one has Eq. (1.3) with Rrad ¼ Jrad ¼ 0 between all other solid materials. The temperature is presumed to be continuous everywhere in the growth apparatus, disregarding the expected step temperature ðkgas grad Tgas Þ . ~ n gas ¼ xb Tgas T ½b at gas– solid interfaces, in absence of precise data for transition coefficients xb : The Stefan–Boltzmann law provides the outer boundary condition (cf. Ref. [14, (3.3)]) 4 ðkgas grad TÞ . ~ Þ; n b ¼ se½b ðTÞðT 4 Troom
ð1:4Þ
s denoting the Boltzmann radiation constant, e denoting the emissivity of the surface, and ~ nb denoting the outer unit normal vector to the solid material b: Condition (1.4) means that the growth apparatus is exposed to a black-body environment
222
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
(e.g. a large isothermal room) radiating at room temperature Troom : Analogously, we use black body phantom closures of open radiation regions such as Gtop and Gbottom depicted in Fig. 1 which emit radiation at Troom ; thereby allowing for radiative interactions between an open cavity and the ambient environment including reflections at the cavity’s surfaces. For the numerical simulations of this work, the above equations will be considered in two space dimensions, taking into account the cylindrical symmetry of the growth system.
2. The induction heating model The axisymmetric crucible is placed inside a copper induction coil as depicted in Fig. 1. An alternating voltage is imposed in the coil, resulting in an alternating current that generates a rapidly oscillating magnetic field, inducing eddy currents in the growth apparatus. The eddy currents cause heat sources due to the Joule effect. We will now describe the modeling of the induction process with the goal of computing the induced eddy currents and the resulting heat sources. It is assumed that a sinusoidal alternating voltage V ðtÞ ¼ V0 sin ðotÞ
ð2:1Þ
is imposed in the coil, since e.g. in experimental 1 setups at the Institut fur (IKZ) in . Kristallzuchtung . Berlin the voltage is more easily measured than the current. To be able to consider the problem in an axisymmetric setting, the actual coil is replaced by N cylindrical rings with a sinusoidal alternating voltage imposed in each ring. Thus, the voltage imposed in the kth ring has the form vk ðtÞ ¼ vk;0 eiot ;
ð2:2Þ
where i denotes the imaginary unit P and v1;0 ; y; vN;0 N are complex voltages satisfying k¼1 vk;0 ¼ V0 : Here, as in the following, we use the complex representation of sinusoidal functions. All solid materials in the growth system are considered as potential conductors, whereas the 1
Institute of Crystal Growth.
gas phase is treated as a perfect insulator. We assume that displacement currents as well as surface currents can be neglected. All the involved electromagnetic quantities are assumed to be axisymmetric and sinusoidal. Using cylindrical coordinates ðr; W; zÞ and their standard basis e W ;~ e z Þ; the complex current density is supposed ð~ e r ;~ to be given by ~ j ¼ jðr; zÞ~ e W : It is shown in Refs. [16,17] that under these assumptions, Maxwell’s equations and Ohm’s law allow one to find a complex-valued magnetic scalar potential fðr; zÞ ~ ¼ eiot fðr; zÞ~ ~ is the magsuch that A e W ; where A netic vector potential, i.e. the magnetic induction ~: ~ ¼ curl A reads B The potential f is determined from the system of elliptic partial differential equations [17, (22), (29)], which we rewrite in the following divergence form (2.3) which is more suitable for our numerical approach via a finite volume discretization: gradðrfÞ n div ¼ 0 in the gas phase; ð2:3aÞ r2 gradðrfÞ i osc f sc vk;0 ¼ þ r2 r 2pr2 in the kth coil ring;
ð2:3bÞ
gradðrfÞ i osc f ¼0 þ r2 r in other conducting materials;
ð2:3cÞ
n div
n div
where n denotes the reciprocal of the magnetic permeability, and sc is the electrical conductivity. The quantities n and sc can vary in space, but they are supposed to be constant in time. System (2.3) is completed by interface and boundary conditions. Owing to the assumption of no surface currents, we have the interface conditions [Ref. 17, Eq. (30)]: npMaterial1 .~ gradðrfÞp n Material1 Material1 r2 npMaterial2 .~ ¼ gradðrfÞp n Material1 Material2 r2 ð2:4Þ on interfaces between Material1 and Material2 ; where p denotes the restriction to the respective material, and ~ n Material1 denotes the outer unit
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
normal vector to Material1 : It is also assumed that f is continuous throughout the whole domain and that f ¼ 0 both on the symmetry axis r ¼ 0 and sufficiently far from the growth apparatus (in Section 3.2, it is found that using a domain radius and height of 1.2 and 1.8 m, respectively, is sufficiently accurate for our purposes). From Ref. [17, Eq. (28)], we get ~ j ¼ j0 eiot~ eW; with j0 ¼ 8 sc vk;0 > < io sc f þ 2pr io sc f > :
in the kth coil ring; in other conducting
ð2:5Þ
materials:
Let Ok denote the two-dimensional domain of the kth coil ring. The total current in the kth ring is then given by eiot Jk ; with Z j0 ðxÞ dx: ð2:6Þ Jk :¼ Ok
Since the coil rings constitute the approximation of a single, connected coil, the total current must be the same in each coil ring, i.e. Jk ¼ Jkþ1
for
k ¼ 1; y; N 1:
ð2:7Þ
Hence, we have to find a complex decomposition V0 ; such that v1;0 ; y; vN;0 of the imposed voltage P Eq. (2.7) is satisfied as well as N v k¼1 k;0 ¼ V0 : This leads to a linear system for v1;0 ; y; vk;0 (see Ref. [18] for more details). After solving the combined system for ðf; v1;0 ; y; vN;0 Þ; we use Eq. (2.5) to compute the power density of the heat sources according to mðt; r; zÞ ¼
jj0 ðr; zÞj2 : 2 sc ðTðt; r; zÞÞ
3. Numerical experiments 3.1. General setup and methods All numerical simulations presented in the following were performed for a growth system consisting of a container with a height of 25 cm and a radius of 8.4 cm placed inside of 5 hollow rectangular-shaped copper induction rings as displayed in Fig. 1, each ring having a surface thickness of 1 mm. The geometric proportions of the coil rings are depicted in Fig. 2. While the horizontal gap between container and coil is fixed at 4.2 cm in all our experiments, the vertical coil position is varied for the simulations in Section 3.3. The prescribed heating voltage satisfies Eq. (2.1). We fix its frequency at f ¼ 10 kHz, whereas we vary the effective voltage such that 190 VpVeff p270 V: The angular frequency and the p amplitude are provided by o ¼ 2pf and V0 ¼ ffiffiffi Veff 2: It is reiterated that only argon is considered in the gas phase and that it is the only material considered a perfect insulator. The solid materials are assumed to be homogeneous and pure materials. In particular, inhomogeneities in the source powder, and chemical changes occurring during the growth process such as sintering and
ð2:8Þ
The assumption of a sinusoidal form of the electromagnetic quantities is only admissible if the involved material parameters are constant in time, and if there is no movement. Thus, to account for the temperature-dependence of the electrical conductivity and for changing coil positions during the transient simulations in Section 3.3, the quasi-stationary electromagnetic problem presented above, including the decomposition of the imposed voltages, is solved in each time step of the transient computations for the heat evolution problem.
223
Fig. 2. Geometric proportions of induction coil rings.
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
224
graphitization of the source, and Si accumulating in the insulation felt are ignored. In real growth systems the coil is usually cooled very efficiently, e.g. by water flowing inside the coil rings. Hence it is reasonable to assume that the coil is always at room temperature, its material data not being changed during the experiments. The complete list of material data used for our numerical simulations of induction heating can be found in Appendix A. The material parameters used for the simulations of the temperature field evolution are identical to those used in Ref. [14] and can be found in Ref. [14, Appendix A]. A finite volume method based on Ref. [19] is used for the space discretizations arising in the stationary computation of the magnetic scalar potential and in the transient temperature simulations, where an implicit Euler scheme provides the time discretization. The discrete scheme has been implemented using the pdelib program package, being developed at the Weierstrass Institute of Applied Analysis and Stochastics (WIAS), Berlin (cf. Ref. [20]). For the numerical methods used for radiation and view factor calculation we refer to Ref. [14, Section 4]. 3.2. Stationary computations of heat sources This section presents and discusses results of three stationary computations at Troom ¼ 293 K of the magnetic scalar potential f and the resulting power densities m: In addition, the total heating ½apparatus powers Ptotal and P½coil total are computed in the growth apparatus and in the coil rings, respectively, by integrating m over the corresponding regions. More precisely Z X ½apparatus ¼ m½b ; Ptotal P½coil total ¼
Z
solids b a coil rings
m½Cu :
solid b
ð3:1Þ
coil rings
In practice the integrals in Eq. (3.1) are approximated by sums over finite volume elements. The simulations use Veff ¼ 230 V, the coil’s upper rim is at z ¼ 14 cm and its lower rim is level with the lower rim of the growth apparatus at z ¼ 0: The results show that even though the
outcome for the magnetic potential distribution does depend on the size of the domain used in the computation, the arising heat source distribution is not affected significantly. The simulations were performed on domains Osmall (radius rsmall ¼ 0:4 m; height hsmall ¼ 1:0 m), Omedium (rmedium ¼ 1:2 m; hmedium ¼ 1:8 m), and Olarge (rlarge ¼ 2:8 m; hlarge ¼ 3:4 m), the growth apparatus being always approximately in the center of the domain. It is noted that instead of imposing the boundary condition f ¼ 0 at r ¼ 0; f ¼ 0 is set at the vertical line r ¼ 1 mm, where f is considered to be sufficiently small. This permits the evaluation of terms with r in the denominator (cf. Eqs. (2.3) and (2.4)). Fig. 3 depicts the results of the three numerical experiments, Osmall in column (a), Omedium in column (b), and Olarge in column (c). Row (1) shows the imaginary part ImðfÞ of the magnetic scalar potential with isolevels spaced at Diso ¼ 105 Wb/m, row (2) shows the real part ReðfÞ (Diso ¼ 107 Wb/m), and row (3) shows the induced power densities m (Diso ¼ 2 MW=m3 ). The pictures demonstrate that ImðfÞ reaches rather far into the surrounding space of the growth apparatus, and the calculated global distribution depends noticeably on the size of the domain. The reach of ReðfÞ is considerably less, as can be seen from the agreement of (2)(b) and (2)(c) in Fig. 3. However, in a close neighborhood of the growth apparatus and the coil rings the distributions of ImðfÞ and ReðfÞ are almost identical for the three domains. This results in virtually identical power density distributions (see row (3) of Fig. 3). Moreover, the values in Table 1 show that the ranges of ImðfÞ; ReðfÞ; and m (both in the coil rings and in the apparatus) agree up to 3 digits between Omedium and Olarge : The same holds true ½apparatus for the total heating powers P½coil : total and Ptotal In consequence Omedium provides a reasonable compromise between accuracy and computing expenditure and is thus used for all computations of heat sources in Section 3.3. The results depicted in Fig. 3 concur qualitatively with the results found in Refs. [9, Fig. 3, 13, Fig. 2]. While a frequency of 10 kHz is not high enough to produce a full-fledged skin effect, the heat
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
225
Fig. 3. Magnetic scalar potential f (1st row: imaginary part, 2nd row: real part) and resulting power density m of heat sources (3rd row) for three different domain sizes. See text for discussion of results.
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
226
Table 1 Extreme values of the real and imaginary parts of the magnetic scalar potential, maximal values of power densities, and values of total powers established during numerical simulations on Osmall ; Omedium ; and Olarge
6
ReðfÞmin ð10 Wb=mÞ ReðfÞmax ð106 Wb=mÞ ImðfÞmin ð104 Wb=mÞ ImðfÞmax ð104 Wb=mÞ 8 3 m½coil max ð10 W=m Þ 7 3 ð10 W=m Þ m½apparatus max P½coil ðkWÞ total P½apparatus ðkWÞ total ½apparatus P½coil þ P ðkWÞ total total
Osmall
Omedium
Olarge
8.80 4.73 1.79 0.00 3.49 1.45 1.63 4.83 6.46
8.82 4.68 1.78 0.00 3.39 1.46 1.57 4.89 6.46
8.82 4.68 1.78 0.00 3.39 1.46 1.57 4.89 6.46
sources decay quickly below a conductor’s surface (see Fig. 3(3), where the maximal power density is established in the lower right-hand corner of the graphite). ½apparatus The values for P½coil in Table 1 total and Ptotal ½apparatus show that some 25 percent of P½coil is total þ Ptotal lost in the coil rings. 3.3. Transient simulation of the temperature evolution In this section we report on transient numerical simulations of the temperature evolution in the growth apparatus, each run starting at T ¼ 293 K. We varied the parameters for argon pressure, heating voltage, and vertical coil positions, to assess their influence during the heating process. It was found in Section 3.2 that the domain Omedium of radius 1.2 m and height 1.8 m is suitable for accurate heat source computations. Hence Omedium is used for all heat source simulations in the current section. As already pointed out at the end of Section 2, the heat source distribution is computed in each time step of the transient simulations to allow for the temperature dependence of electrical conductivities (cf. Eqs. (A.1)). The temperature distribution is only calculated in the growth apparatus, i.e. on a much smaller grid GT than the grid Gm used for the heat source computations. In consequence in each time step data for T are carried from GT to Gm and data for m are carried from Gm to GT :
It is noted that the relaxation times of the pressure and of the temperature in the gas phase always lie orders of magnitude below the relaxation times of the temperature in the solid components. Since the gas does not have a significant influence on the temperature of the solid parts, each temporal snapshot of a transient simulation furnishes the quasi-stationary temperature distribution in the gas phase determined by the temperature distribution on the adjacent solid surfaces. In particular, for a constant pressure pðArÞ ; the time derivative in Eq. (1.1b) vanishes identically as can be seen from Eq. (1.2) and the gas law pðArÞ ¼ rðArÞ ðR=M ðArÞ ÞT: Thus, the quasistationary temperature distribution in the gas phase does not depend on the gas pressure. The argument of the preceding paragraph suggests that Eq. (1.1a) could be coupled with the stationary version of Eq. (1.1b), i.e. with divðkðArÞ grad TÞ ¼ 0: However, for the numerical regularity of the problem, it is advantageous to keep the time derivative in the gas equation which is also helpful within the present code implementation. For the numerical experiments described in the following, instead of using a given argon pressure, we fixed rðArÞ ¼ 3:73 103 kg=m3 ; which corresponds to T ¼ 2575 K and pðArÞ ¼ 2 103 Pa. The total heating power in the growth apparatus P½apparatus and in the coil rings P½coil total total changes during each simulation run due to the temperature-dependent electrical conductivities. Moreover, if present, a movement of the induction coil also changes P½apparatus and P½coil total total : For one of our simulations, using Veff ¼ 230 V and the fixed coil position C16 2 according to the notation introduced in the next paragraph, Figs. 4(a)(1) and (a)(2) depict the time evolution of P½apparatus and P½coil total total , respectively. The curves show an inverse behavior, i.e. P½apparatus decreases as P½coil total total increases and vice versa. The peaks after some 500 s in Fig. 4(a) correspond to the peaks in Fig. 4(b) which show the electrical conductivity of the graphite crucible depending on temperature in Fig. 4(b)(1) and depending on the maximal temperature established in the graphite at time t in Fig. 4(b)(2). The minimum in Fig. 4(a)(1) matches the maximum in Fig. 4(b)(2), showing that the system is in a state
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
227
Fig. 4. Time evolution of total heating power in growth apparatus (a)(1) and in coil rings (a)(2) for numerical experiment C16 2 230 V in comparison with the temperature-dependent electrical conductivity of the graphite crucible (b).
where increasing the electrical conductivity of the apparatus (lowering its resistance) reduces ½apparatus Ptotal : We now discuss three series of numerical 16 18 experiments, referred to as C14 0 ; C2 ; and C4 ; respectively, each series using a different vertical coil position. For C14 the coil is positioned 0 between z ¼ 0 and 14 cm, i.e. the lower rim of the bottom coil ring is at z ¼ 0; and the upper rim of the top coil ring is at z ¼ 14 cm. The meaning of 18 C16 2 and C4 is analogous. There are five experiments in each series using different values for Veff ; namely Veff ¼ 190, 210, 230, 250, and 270 V. Subsequently, a fourth series C18-14 4-0 is considered, using a moving induction coil. For each simulation, Table 2 contains P½apparatus and P½coil total total ; i.e. the
total heating power in the growth apparatus and in the coil rings, respectively, established in the stationary final state. Results for the time evolution and the final values of the temperature evolution at four points of particular importance are depicted in Figs. 5–7, respectively. The monitored points are located at the blind holes at the top and at the bottom of the growth apparatus, and at the surface of the crystal seed and of the source powder, respectively (labeled Ttop ; Tbottom ; Tseed ; and Tsource in Fig. 1). The significance of Ttop and Tbottom lies in the blind holes being the principal locations for temperature measurements; Tseed and Tsource are of interest, since their difference is a key control factor for the crystal’s growth rate and quality (cf. e.g.
228
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
Fig. 5. Time evolution of Ttop (row (1)), Tbottom (row (2)), Tseed (row (3)), and Tsource Tseed (row (4)) for the simulation series C14 0 (column (a)) and C16 2 (column (b)).
Refs. [2,5]). Results of the time evolution of the global temperature field and the heat sources are shown in Fig. 8, and results of the quasi-stationary
temperature distributions at the crystal and its immediate surroundings are depicted in Fig. 9. We now proceed to discuss our results in detail.
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
229
Fig. 6. Time evolution of Ttop (row (1)), Tbottom (row (2)), Tseed (row (3)), and Tsource Tseed (row (4)) for simulation series C18 4 (column 18-14 (a)) and C4-10 (column (b)).
We find that for each time t; the temperatures Ttop ; Tbottom ; Tseed ; and Tsource depend nearly linearly on the heating voltage Veff : For Ttop ;
Tbottom ; and Tseed ; this can be observed in Figs. 5 and 6. Moreover, for the quasi-stationary state at t ¼ 30 000 s, the nearly linear dependence of these
230
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
Fig. 7. Temperatures Ttop (in (1)), Tbottom (in (2)), Tseed (in (3)), and Tsource Tseed (in (4)) for the quasi-stationary state at t ¼ 30 000 s 16 18 for the simulation series C14 0 ; C2 ; and C4 depending on the effective heating voltage.
temperatures and also of Tsource Tseed is illustrated in Fig. 7. Figs. 5–7 show furthermore that Ttop ; Tseed ; and Tsource increase significantly be18 tween C14 0 (low coil position) and C4 (high coil position), whereas Tbottom remains almost constant. This can be explained by the heat sources being shifted upwards and the source powder’s insulating property which, though less prominent, is still effective even at high temperatures (cf. Fig. 8(a) and Ref. [14, (A.6b)]), thereby hindering heat generated below the source powder from reaching regions above the source powder.
The shape of the curves in row (4) of Figs. 5 and 6 depicting Tsource Tseed is also due to the insulating characteristic of the source powder and the gas phase, initially keeping the powder’s surface below the temperature of the seed which the heat reaches via conduction through the graphite. This causes the negative peak in the corresponding graphs. At higher temperatures, radiative heat transfer becomes more effective both in the gas phase and the porous source, resulting in the source growing warmer than the seed, as is required for crystal growth. Moreover,
18-14 Fig. 8. Time evolution of heating process for numerical experiment C4-0 230 V (coil moving at 1:33 cm/h). Column (a): temperature evolution, difference between neighboring isotherms is 20 K. Column (b): heat source evolution, difference between isolevels is 10 kW/m3 ; darker regions indicate larger heat sources.
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
231
232
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
Fig. 9. Close-ups of the temperature distribution in the crystal, the gas phase, and the crucible directly above the crystal for the quasistationary states at t ¼ 30 000 s of the numerical experiments with lowest and medium coil position in (a) and (b), respectively. Temperature difference between neighboring isotherms is 5 K.
it is seen in Fig. 7(4) that at the quasi-stationary state, a lower coil position or a higher heating voltage result in Tsource Tseed being increased as heat reaches the powder’s surface more readily at higher temperatures due to radiative heat transfer. Comparing our findings in Fig. 7(4) with [21, Fig. 6] establishes a qualitative agreement, whereas a quantitative comparison is meaningless due to the different setups. Fig. 7(4) indicates that Tsource Tseed can be tuned effectively by controling the heating voltage and the vertical coil position. We now discuss the series of numerical experiments C18-14 4-0 ; where the coil starts at the position used for C18 4 ; moves downwards at the rate 1.33 cm/h, and stops at the position used for C14 0 : During these simulations the grid for the heat source computations is newly generated in each time step. The objective of the series C18-14 is to 4-0 study how a moving coil affects the temperature field evolution during the heating process. In this article, we do not investigate coil movements that follow the growing crystal in the quasi-stationary temperature field which are often used in growth experiments to compensate the influence of the growing crystal on the temperature distribution. The rate of such coil movements is usually in the order of 1 mm/h. As expected, the curves for the moving coil (Fig. 6(b)) initially coincide with the corresponding curves of C18 4 (Fig. 6(a)) and coincide with the corresponding curves of C14 0 (Fig. 5(a)) after the coil’s lowest position is reached at t ¼ 3 h: For experiment C18-14 4-0 230 V, we discuss our results concerning the time evolution of the global temperature distribution and the heat sources. It can be seen in column (a) of Fig. 8 that the minimal temperature Tmin is always established at
the outside of the outer insulation material, where the isotherms become very dense at higher temperatures, producing the dark outer regions in Fig. 8(2)(a) and (3)(a). The maximal temperature Tmax is found in the graphite strip between source and insulation in Fig. 8(1)(a) and inside the labeled isotherms below the source in Fig. 8(2)(a) and (3)(a). The change of the location of Tmax corresponds to the heat sources moving downwards together with the induction coil according to Fig. 8(b). Fig. 8(a) also demonstrates how the insulating property of the gas phase and the source powder changes with time, causing a large temperature gradient in the gas phase in Fig. 8(1)(a) and local minima in the source in Figs. 8(1)(a) and (2)(a). Due to radiative heat transfer through the gas phase and through the pores of the source, the temperature on the boundary of the growth chamber (and hence also inside the gas itself) is more homogeneous in Figs. 8(2)(a) and (3)(a). Moreover, the local minimum has vanished in Fig. 8(3)(a), even though some of the source’s insulating quality persists. The increase in P½apparatus between Figs. 8(1)(b) total and (2)(b) followed by the drop between Figs. 8(2)(b) and (3)(b) is consistent with Fig. 4(a)(1) and Table 2: The increase corresponds to the rise after the initial minimum in Fig. 4(a)(1), and the drop matches the power’s decrease with a lower coil position found in Table 2. Finally, Fig. 9 displays close-ups of the final temperature distributions in the crystal’s immediate surroundings for two coil positions 16 C14 0 230 V (in Fig. 9(a)) and C2 230 V (in Fig. 9(b)). The picture in Fig. 9(a) also constitutes a close-up of Fig. 8(3)(a), since the stationary states are identical for the simulations C14 0 230 V
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
233
Table 2 Total heating power P½apparatus in the growth apparatus and P½coil total total in the coil rings for the quasi-stationary state at t ¼ 30 000 s, depending on heating voltage and coil positions Series
Veff P½apparatus ðkWÞ total
C14 0 C16 2 C18 4 C18-14 4-0
P½coil total ðkWÞ
190
210
230
250
270
190
210
230
250
270
3.15 3.40 3.55 3.15
3.87 4.19 4.37 3.87
4.67 5.06 5.27 4.67
5.55 6.01 6.26 5.55
6.51 7.05 7.34 6.51
1.08 1.09 1.10 1.08
1.32 1.33 1.37 1.32
1.58 1.59 1.60 1.58
1.86 1.88 1.89 1.86
2.17 2.18 2.20 2.17
and C18-14 4-0 230 V. The lower temperatures found for the lower coil position in Fig. 9, are once more explained by the powder source forming a barrier for heat generated in the lower part of the apparatus. Moreover, Fig. 9 shows that the final temperature gradient between source and seed is nearly independent of the coil position and approximately constant along the r ¼ 0 axis with a slight increase at the seed.
. Bottcher, Thomas Muller, . Detlev Schulz, and Dietmar Siche. We also thank Jurgen . Sprekels and Manfred Uhle of the WIAS in Berlin for helpful discussions and advice. This research was partly funded by the Bundesministerium fur . Bildung, Wissenschaft, Forschung und Technologie2 (BMBF) within the programs Neue Mathematische Verfahren in Industrie und Dienstleistungen3 # 03SPM3B5 and Mathematische Verfahren zur Losung von Problemstellungen in . Industrie und Wirtschaft 4 # 03SP7FV16.
4. Conclusions The influence of induction heating during PVT was studied by means of transient numerical simulations including diffusive and radiative heat transfer, where a quasi-stationary elliptic problem is solved in each time step to determine temperature- and coil-position-dependent heat sources. We presented results of computed temperature evolutions depending on heating voltage and vertical coil position. The results agree with physical expectations and provide insights on how to control both the global temperature distribution evolution and the temperature evolution at the seed and the source by adjusting heating voltage and vertical coil position.
Appendix A. Material data for induction heating
Acknowledgements
2 German Ministry for Education, Science, Research, and Technology. 3 New Mathematical Methods in Manufacturing- and Service Industry. 4 Mathematical Methods Solving Problems in Industry and Business.
The IKZ in Berlin has provided us with material data and experimental data throughout the project. In particular, we wish to thank Klaus
The relevant materials are argon, which is the only material treated as a perfect insulator, copper, graphite crucible, insulation, SiC powder, and SiC crystal. In each material the reciprocal of the magnetic permeability is assumed to be n ¼ m1 0 ¼ ð107 =4pÞ A m=V s: A.1. Mass densities of conducting materials The values of Table 3 are according to Refs. [22 (Cu); 23 (graphite crucible, insulation); 24 (SiC powder); 25 (SiC crystal)].
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235
234
Table 3 Mass densities of conducting materials Material
Cu
Graphite crucible
Insulation
SiC powder
SiC crystal
r ½kg=m3
8930
1750
170
1700
3140
A.2. Electrical conductivities of conducting materials According to Ref. [22] the value for the electrical conductivity of copper at room temperature is s½Cu ¼ ð1:7 108 O mÞ1 ¼ 5:9 107 ð1=O mÞ: The material parameters used for the electrical conductivities of the remaining conducting materials are s½Crucible ðTÞ ¼ !!1 6 ln ðT=1023 KÞ 2 10 ; 28:9 18:8 exp 2:37 Om ðA:1aÞ s½Insulation ðTÞ ¼
245:7 1 ; 1 þ T=2500 K O m
s½SiCPowder ðTÞ ¼ 100
1 ; Om
ðA:1bÞ ðA:1cÞ
1 : ðA:1dÞ Om Eq. (A.1a) is according to Ref. [26], and Eqs. (A.1b), (A.1c), and (A.1d) are according to Ref. [27]. s½SiCCrystal ðTÞ ¼ 105
References [1] St.G. Muller, . R.C. Glass, H.M. Hobgood, V.F. Tsvetkov, M. Brady, D. Henshall, J.R. Jenny, D. Malta, C.H. Carter Jr., J. Crystal Growth 211 (2000) 325. [2] A.O. Konstantinov, Sublimation growth of SiC, in: G.L. Harris (Ed.), Properties of Silicon Carbide, EMIS Datareview Series, no. 13, Institution of Electrical Engineers, INSPEC, London, pp. 170–203. [3] D.L. Barrett, J.P. McHugh, H.M. Hobgood, R.H. Hopkins, P.G. McMullin, R.C. Clarke, W.J. Choyke, J. Crystal Growth 128 (1993) 358. . [4] H.-J. Rost, D. Siche, J. Dolle, W. Eiserbeck, T. Muller, D. Schulz, G. Wagner, J. Wollweber, Mater. Sci. Eng. B 61–62 (1999) 68.
[5] N. Schulze, D.L. Barrett, G. Pensl, Appl. Phys. Lett. 72 (13) (1998) 1632. [6] A.S. Segal, A.N. Vorob’ev, S.Yu. Karpov, E.N. Mokhov, M.G. Ramm, M.S. Ramm, A.D. Roenkov, Yu.A. Vodakov, Yu.N. Makarov, J. Crystal Growth 208 (2000) 431. [7] S.Yu. Karpov, Yu.N. Makarov, M.G. Ramm, Phys. Stat. Sol. B 202 (1997) 201. [8] K. Chourou, M. Anikin, J.M. Bluet, J.M. Dedulle, R. Madar, M. Pons, E. Blanquet, C. Bernard, P. Grosse, C. Faure, G. Basset, Y. Grange, Mater. Sci. Eng. B 61–62 (1999) 82. [9] M. Pons, M. Anikin, K. Chourou, J.M. Dedulle, R. Madar, E. Blanquet, A. Pisch, C. Bernard, P. Grosse, C. Faure, G. Basset, Y. Grange, Mater. Sci. Eng. B 61–62 (1999) 18. [10] 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, Yu.N. Makarov, Mater. Sci. Eng. B 61–62 (1999) 107. [11] S.Yu. Karpov, A.V. Kulik, I.A. Zhmakin, Yu.N. Makarov, E.N. Mokhov, M.G. Ramm, M.S. Ramm, A.D. Roenkov, Yu.A. Vodakov, J. Crystal Growth 211 (2000) 347. [12] M. Selder, L. Kadinski, Yu. Makarov, F. Durst, P. Wellmann, T. Straubinger, D. Hoffmann, S. Karpov, M. Ramm, J. Crystal Growth 211 (2000) 333. [13] Q.-S. Chen, H. Zhang, V. Prasad, C.M. Balkas, N.K. Yushin, A System Model for Silicon Carbide Crystal Growth by Physical Vapor Transport Method, 1999 National Heat Transfer Conference, NHTC99-222, ASME, 1999, pp. 1–8. ! [14] O. Klein, P. Philip, J. Sprekels, K. Wilmanski, J. Crystal Growth 222 (2001) 832. ! [15] N. Bubner, O. Klein, P. Philip, J. Sprekels, K. Wilmanski, J. Crystal Growth 205 (1999) 294. [16] S. Clain, J. Rappaz, M. Swierkosz, Coupling between nonlinear Maxwell and heat equations for an induction heating problem: modeling and numerical methods, in: M. Krizek et al. (Ed.), Finite element methods. 50 years of the Courant element, Lecture Notes in Pure and Applied Mathematics Vol. 164, Marcel Dekker, New York, NY, pp. 163–171. [17] J. Rappaz, M. Swierkosz, Modelling in numerical simulation of electromagnetic heating, in: Modelling and Optimization of Distributed Parameter Systems, Warsaw, 1995, Chapman & Hall, New York, pp. 313–320. [18] O. Klein, P. Philip, IEEE Trans. Mag. 38 (3) (2002) 1519. [19] J. Fuhrmann, On numerical solution methods for nonlinear parabolic problems, in: Modeling and Computation
O. Klein, P. Philip / Journal of Crystal Growth 247 (2003) 219–235 in Environmental Sciences, Notes on Numerical Fluid Mechanics, Vol. 59, Vieweg, Braunschweig, Germany, pp. 170–180. [20] J. Fuhrmann, Th. Koprucki, H. Langmach, pdelib: an open modular tool box for the numerical solution of partial differential equations. Design patterns, Proceedings of the 14th GAMM Seminar on Concepts of Numerical Software, Kiel, January 23–25, 1998 (Kiel, Germany), 2001. [21] D.L. Barrett, R.G. Seidensticker, W. Gaida, R.H. Hopkins, W.J. Choyke, J. Crystal Growth 109 (1991) 17. [22] G. Benkowsky, Induktionserw.armung. H.arten, Gluhen, . . Schmelzen, Loten, Schweien, 5th edition, Verlag Technik GmbH, Berlin, 1990.
235
[23] T. Muller, . D. Schulz, D. Siche, Material vendor’s data, communicated by the Institut fur . Kristallzuchtung, . Berlin, 1999. [24] T. Muller, . D. Schulz, D. Siche, Institut fur . Kristallzuchtung, . Berlin, 1999, unpublished experimental data. [25] O. Nilsson, H. Mehling, R. Horn, J. Fricke, R. Hofmann, S.G. Muller, . R. Eckstein, D. Hofmann, High Temp. – High Press. 29 (1) (1997) 73. [26] J.G. Hust (Ed.), Standard Reference materials: finegrained, isotropic graphite for use as NBS thermophysical property RM’s from 5 to 2500 K, Tech. Report NBS/SP260/89, National Bureau of Standards, September 1984. [27] P. R(aback, Modeling of the sublimation growth of silicon carbide crystals, Ph.D. Thesis, Helsinki University of Technology, 1996.