- Email: [email protected]
Global simulation of a silicon Czochralski furnace
Journal of Crystal Growth 234 (2002) 32–46
Global simulation of a silicon Czochralski furnace Mingwei Lia, Yourong Lia, Nobuyuki Imaishia,*, Takao Tsukadab a
Institute of Advanced Material Study, Kyushu University, 6-1 Kasuga-koen, Kasuga, Fukuoka 816-8580, Japan Institute of Chemical Reaction Science, Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai 980-8577, Japan
b
Received 6 July 2000; accepted 25 July 2001 Communicated by J.J. Derby
Abstract To understand the characteristics of the Czochralski (Cz) furnace for the single-crystal growth of silicon, a set of global analyses of momentum, heat and mass transfer in small Cz furnaces (crucible diameter: 7.2 cm, crystal diameter: 3.5 cm, operated in a 10 Torr argon flow environment) is carried out using the finite-element method. The global analysis assumes a pseudosteady axisymmetric state with laminar flow, equilibrium relations at the melt/silica interface and vapor–liquid chemical equilibrium at the melt/gas interface, as well as the segregation coefficient of unity at the melt/crystal interface. Convective and conductive heat transfers, radiative heat transfer between diffuse surfaces and the Navier–Stokes equations for gas and melt phases are all combined and solved together. Thus, the velocities and temperatures obtained are used to calculate the oxygen concentrations. The global analysis code is effectively used to discuss the influences of the Marangoni effect and a gas guide (or a heat shield) placed between the crucible and the crystal. The results indicate that the gas guide reduces the heater power and changes the melt flow pattern and oxygen transport. The melt flow pattern is strongly dependent on the Marangoni effect and gas flow near the surface, and changes the oxygen concentration significantly. This analysis reveals the importance and effectiveness of global analysis. r 2002 Elsevier Science B.V. All rights reserved. Keywords: A1. Computer simulation; A1. Fluid flows; A1. Heat transfer; A1. Mass transfer; A2. Czochralski method; B2. Semiconducting silicon
1. Introduction For the production of high-quality silicon crystals by the Czochralski (Cz) method, the control of concentration and uniformity of oxygen, as well as of the temperature gradient, in the growing crystal is required. In order to control these factors, most of the recent Cz furnaces are equipped with substructures in the furnace, such as *Corresponding author. Tel: +81-92-583-7793; fax: +81-92583-7796. E-mail address: [email protected] (N. Imaishi).
a funnel-shaped radiation shield and multiple heaters for the sake of thermal environment control as well as for saving energy. To control oxygen concentration, the operating pressure and the gas flow rate are included as operating parameters, in addition to conventional control factors such as crystal and crucible rotation rates. The application of various types of magnetic fields to the electrically conducting silicon melt is fast becoming popular for stabilizing turbulent melt flow and controlling oxygen transfer in large-scale crucibles. The oxygen in the crystal originates from the dissolution of silica from the crucible into
0022-0248/02/$ - see front matter r 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 0 1 ) 0 1 6 3 4 - 7
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
molten silicon. Dissolved oxygen is convectively transported to the melt/crystal interface and incorporated into the crystal. It is believed that before arriving at the melt/crystal interface, most of the dissolved oxygen is evaporated into the gas phase through the melt/gas interface via the chemical reaction SiðlÞ þ OðlÞ"SiOðgÞ at the melt/gas interface. In recent years, many numerical analyses [1–7] have been conducted on melt flow and oxygen transport in the silicon melt. These analyses have helped us greatly in understanding the transport mechanism. The experimental results obtained by Kawanishi et al. [8] and Kanda et al. [9] suggest that the melt surface flow plays a very important role in the radial nonuniformity of oxygen concentration in crystal. However, the boundary condition for the diffusion equation at the melt/gas interface in these numerical simulations assumed either a zero oxygen concentration [1,3,4] or a constant evaporation coefficient [5–7]. The zero oxygen concentration boundary condition implies that there is no diffusion resistance in the gas phase. This is in contrast to the experimental results obtained by Machida et al. [10], which revealed that the gas phase pressure and the gas flow rate significantly affect the oxygen concentration in the crystal. In their experiments, the increase of gas pressure or decrease of gas flow rate caused a considerable increase of oxygen concentration in the crystal. In a crystal grown in a Cz furnace with a heat shield, which also acts as a gas guide, the oxygen concentration decreases until a critical gas flow rate is reached. Beyond the critical flow rate, however, the oxygen concentration increases again with further increase of the gas flow rate. Machida et al. speculated that at higher gas flow rates, the shear stress from the gas flow may change the melt flow pattern, decrease the evaporation rate of SiO and consequently increase the oxygen concentration in the grown crystal. Bornside et al. [11] conducted the first numerical simulation of the evaporation of SiO and the incorporation of carbon through the gas phase in a conventional Cz furnace. However, in their analysis, they assumed that the melt surface was rigid with a constant vapor pressure of SiO. The effect of substructures in the Cz furnace, such as a heat shield, has not yet been analyzed,
33
except for some analyses on their effect on the melt/crystal interface shape [12,13]. The basic structure of the Cz furnace is very simple and the number of components is rather small. Despite such a simple structure, phenomena in the Cz furnace are highly nonlinear and they are strongly coupled with each other. Therefore, one small change may cause various effects. For example, a cylindrical or funnel-shaped heat shield changes the temperature gradient in the crystal as well as the shape of the melt/crystal interface [12,13]. At the same time, the heat shield changes the path of gas flow near the melt/gas interface as well as the evaporation rate of oxygen. The heat shield reduces the radiative heat loss and in turn, the heater power, temperature distributions and melt flow are also modified. These change the solubility of SiO2 in silicon melt, the equilibrium vapor pressure of SiO at the melt/gas interface, as well as the gas-phase mass-transfer coefficient. These are all closely related to each other and no individual component analysis can claim its own significance. In order to understand the overall characteristics of the Cz furnace, ‘‘global analysis’’ must be applied. Previously, Kinney and Brown [4] and Muhe et al. [14] analyzed the melt convection, and conductive and radiative heat transfer in a Cz furnace, but did not pay any attention to the gasphase transport phenomena. In the present work, a series of simulations was conducted for a small silicon Cz furnace using a laboratory-made global analysis code that includes heat exchange in the entire furnace, melt and gas flows, and oxygen and SiO transport in the melt and in the gas. The code requires input data for only one set of operating conditions, such as gas flow rate, gas-phase pressure, heater power, furnace chamber wall temperature, geometry of the furnace, diameter and length of the crystal, and melt volume. The pseudosteady state and axisymmetry were assumed. The melt/crystal interface shape was determined as the isotherm at the melting temperature that satisfies the heat balance as well as all boundary conditions. The effects of a gas guide with a simple shape were discussed, together with the Marangoni (thermocapillary) effect. For simplicity, crystal and crucible rotations were neglected in this work.
34
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
2. Model formulation 2.1. Basic assumptions The following assumptions are introduced to the model: (1) The flow is laminar in both the melt and the gas and the melt is incompressible. (2) Thermophysical properties of the melt are assumed to be constant except for the temperature dependence of density and surface tension, in the calculation of the buoyancy force and the Marangoni effect, respectively. However, the thermophysical properties of argon gas are treated as temperature-dependent. (3) The PVT relation of argon gas is expressed by the ideal gas law. (4) No supercooling occurs at the melt/crystal interface and the interface shape coincides with the melting temperature isotherm. (5) Oxygen concentration in the argon gas at the inlet is zero. (6) Oxygen concentration at the melt/crucible interface is given by the solubility equilibrium [15]. (7) The oxygen concentration and the partial pressure of SiO are in chemical equilibrium at the melt–gas interface [11]. (8) The partial pressure of oxygen in the bulk of argon gas is zero. (9) The segregation coefficient of oxygen is assumed to be unity. (10) The temperature coefficient of surface tension is taken as a disposable parameter to test the role of the Marangoni effect. (11) Evaporated SiO does not deposit on any solid surface. Among these assumptions, only the last one is not realistic but may be acceptable as long as the discussion is limited to the evaporation rate of SiO. For further analysis, the vapor pressure of SiO and its reaction with carbon materials should be included, as in Ref. [11]. Assumption (9) is introduced here because there is no reliable chemical analysis data on the argon gas used in the experiments. If such data were available, the oxygen supply from the gas phase must be taken into account, as in Ref. [16]. 2.2. Mathematical models 2.2.1. Governing equations Under the above assumptions, the general governing equations for flow, temperature and concentration fields are given as follows.
In the melt: r vl ¼ 0;
ð1Þ
rl vl rvl ¼ rpl r sl þ rl gbðTl Tm Þez ;
ð2Þ
rl Cpl vl rTl ¼ ll r rTl ;
ð3Þ
vl rxO ¼ DO r rxO :
ð4Þ
In the gas: r ðrg vg Þ ¼ 0;
ð5Þ
rg vg rvg ¼ rpg r sg rg gez ;
ð6Þ
rg Cpg vg rTg ¼ r ðkg rTg Þ;
ð7Þ
r ðcg vg ySiO Þ ¼ r ðcg DSiO rySiO Þ:
ð8Þ
In the crystal: rs Cps Vs ez rTs ¼ r ðls rTs Þ:
ð9Þ
In the heater: r ðlh rTh Þ þ Q ¼ 0:
ð10Þ
In the other solid materials: r ðli rTi Þ ¼ 0:
ð11Þ
T; v; p; and s are the temperature, velocity vector, pressure and stress tensor, respectively. Tm is the melting temperature, r the density, l the thermal conductivity, D the mass diffusivity, b the thermal expansion coefficient, g the gravitational acceleration constant, Cp the heat capacity, Vs the crystal pulling rate and c the molar density. x and y are the mole fractions of oxygen and SiO, respectively. Subscripts l; g; s; h; i; O and SiO indicate melt, gas, crystal, heater, other materials, oxygen and SiO, respectively. ez is a unit vector parallel to the z axis. 2.2.2. Boundary conditions for the temperature field At the melt/gas interface: ll n rTl ¼ lg n rTg þ qrad;l :
ð12aÞ
At the melt/crystal interface: ll n rTl ls n rTs ¼ rs Vs DHs n ez ;
ð12bÞ
Tl ¼ Ts ¼ Tm :
ð12cÞ
35
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
At the crystal side and top: ls n rTs ¼ lg n rTg þ qrad;s : At the other exposed surfaces: li n rTi ¼ lg n rTg þ qrad;i :
ð12eÞ
At the inlet of Ar gas: Tg ¼ Tg0 :
ð12fÞ
At the outlet of gas: n rTg ¼ 0:
tl : ney tg : ney ¼ gT rTl ey ;
ð14kÞ
ey vl ey vg ¼ 0:
ð14lÞ
ð12dÞ
ð12gÞ
In these equations, qrad;i is the net radiative heat flux per unit area on the surfaces, and is given by " ! # X 4 4 qrad;i ¼ ei sB Ti Aj ej Tj Gji =Ai =ei ; ð13Þ
At the melt/crystal interface: u1 ¼ 0;
ð14mÞ
v1 ¼ 0;
ð14nÞ
wl ¼ rOs :
ð14oÞ
At the other surfaces of crystal and the puller surface: ug ¼ 0;
ð14pÞ
vg ¼ 0;
ð14qÞ
wg ¼ rOs :
ð14rÞ
j
where Gji is Gebhart’s absorption factor, which is the fraction of emission from surface Aj that reaches Ai and is absorbed. Gji as well as the view factors between any two surface elements in the furnace is calculated using our code [13]. DHs is the latent heat of fusion, n the unit normal vector of the interface element and e the emissivity. 2.2.3. Boundary conditions for the flow field At the side wall and the bottom of the crucible:
At the inlet of Ar gas: ug ¼ 0;
ð14sÞ
wg ¼ 0;
ð14tÞ
vg ¼ Vin :
ð14uÞ
At the outlet of Ar gas: t n ¼ 0:
ð14vÞ
u1 ¼ 0;
ð14aÞ
v1 ¼ 0;
ð14bÞ
ug ¼ 0;
ð14wÞ
wl ¼ rOc :
ð14cÞ
vg ¼ 0;
ð14xÞ
At the other crucible surfaces and the pedestal surface:
wg ¼ 0:
ð14yÞ
ug ¼ 0;
ð14dÞ
vg ¼ 0;
ð14eÞ
wg ¼ rOc :
ð14fÞ
u; v; and w are the radial, axial and the azimuthal components of velocity, and Os and Oc are the angular rotation rates of the crystal and the crucible, respectively. t and ey are the unit tangent and azimuthal vectors of each interface element, respectively. The melt/crystal interface shape is determined such that Eq. (12c) is satisfied, i.e., the interface coincides with the melting-point isotherm.
At the melt/gas interface:
At other motionless solid surfaces:
n vl ¼ 0;
ð14gÞ
n vg ¼ 0;
ð14hÞ
tl : nt tg : nt ¼ gT rTl t;
ð14iÞ
2.2.4. Boundary conditions for oxygen and SiO At the crucible wall [15]:
t vl t vg ¼ 0;
ð14jÞ
xO ¼ 0:01eð7:822 370=TÞ :
ð15aÞ
36
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
At the melt/gas interface [11]: ySiO ¼ 101 325xO e
ð17:821 000:0=TÞ
=P0 ;
cg DSiO n rySiO ¼ cl DO n rxO :
2.4. Simulation conditions ð15bÞ ð15cÞ
At the melt/crystal interface: n rxO ¼ 0:
ð15dÞ
At the inlet of Ar gas: ySiO ¼ 0:
ð15eÞ
At the outlet of gas: qySiO ¼ 0: qz At the other exposed surfaces:
ð15fÞ
n rySiO ¼ 0;
ð15gÞ
where P0 (=1333 Pa or 10 Torr) is the total pressure. 2.3. Numerical procedure To solve this problem, the Galerkin finiteelement method is used. The calculation domain is discretized by isoparametric quadrilateral elements. In each element, velocity vectors, temperature and concentration are approximated with the bilinear trial function f as follows: X vðr; zÞ ¼ f i vi ; ð16aÞ X ð16bÞ Tðr; zÞ ¼ fi Ti ; X ð16cÞ yðr; zÞ ¼ fi ySiOi ; X ð16dÞ xðr; zÞ ¼ fi xOi : The pressure is assumed to be constant in each element. By substituting Eqs. (16a)–(16d) into the governing equations and the boundary conditions for the flow, temperature and concentration fields, the Galerkin procedure gives a set of algebraic equations. The melt/gas interface shape is calculated by solving the Young–Laplace equation using the Galerkin method with the Hermit cubic trial function. In the following analyses, this melt/gas interface shape is approximated by 50 linear segments.
The configurations of the Cz furnaces adopted in the present global analyses are shown in Fig. 1. The inner diameter of both crucibles is 7.2 cm, and the diameter of both crystals is 3.5 cm. The heat shield (or gas guide) is a straight cylinder (inner diameter 5.15 cm, thickness 2.0 mm) and its lower end is located 1.14 cm above the melt surface, as shown in Fig. 1. It is unrealistic to install a heat shield in such a small furnace. The heat shield is adopted in this analysis to elucidate the distinguishable effects of such a substructure on heat transfer and mass transfer, and also to demonstrate the importance of the global analysis. In this paper, crystal and crucible rotations, as well as any azimuthal distribution, are neglected. Thus, we assume Os ¼ Oc ¼ 0 and neglect Eqs. (14k) and (14l). In the present work, the diameter and length of the crystal and the heater power are given. The crystal pull rate is treated as one of the unknown
Fig. 1. Configuration of the Cz furnaces used in this study.
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
variables, together with velocities, pressures, temperatures, and the melt/crystal interface coordinates, by simultaneously solving the set of nonlinear algebraic equations using the Newton– Raphson method. The total number of elements, nodal points and unknowns are 54 542, 55 107 and 242 486, respectively. The concentration field is computed by using much finer elements, especially in the melt phase, with 81 576 elements, 82 458 nodal points and 82 458 unknown variables. The sizes of the minimum elements that are located at the corner of the melt/gas interface and the crucible wall in the melt are Drmin ¼ 0:125 mm, Dzmin ¼ 0:073 mm for the flow and heat transfer analysis, and Drmin ¼ 2:60 mm, Dzmin ¼ 0:71 mm for the oxygen transport analysis. For mass transfer calculations, the temperature and velocities at each point are interpolated based on the obtained velocity and temperature values. The code requires only a set of operating conditions as input data, such as gas flow rate, gas phase pressure, heater power, furnace chamber wall temperature, geometry of the furnace, diameter and length of the crystal, and total volume of silicon (melt and crystal). Typical conditions are listed in Table 1, together with the thermophysical properties used in this work. The temperature coefficient of surface tension is assumed to be either zero or 7 105 N/mK. The latter value is taken from the recent experiments of Nogi [17] (In their final report [18], the value is revised to – 6.2 105 N/mK.). The heater power distribution is evaluated by solving the Laplace equation for the current in the heater, and then the power distribution is averaged over both azimuthal and radial directions expressed as a function of the longitudinal coordinate.
37
Table 1 Physical properties and processing parameters Thermal conductivity l ¼ 22:0 W/mK (crystal), l ¼ 64U0 W/mK (melt), l ¼ 2:98 W/mK (quartz) l ¼ 60:0 W/mK (graphite, heater, pedestal, radiation shield and gas guide) l ¼ 59:0 W/mK (puller) Viscosity m ¼ 7U0 104 kg/ms (melt) Diffusivity DO ¼ 5 108 m2/s, DSiO ¼ 8:626 106 T 1:75 =P0 m2/s Density r ¼ 2530:0 kg/m3 (melt), r ¼ P0 M=RT kg/m3 (Argon gas) Thermal expansion coefficient of melt bl ¼ 1:5 104 K1 Surface tension and its temperature coefficient g ¼ 0:735 N/m, gT ¼ 0:0 or 7.0 105 N/mK Heat of fusion DHm ¼ 1410 J/kg Heat capacity Cpl ¼ 1000 J/kgK, Cps ¼ 1000 J/kgK, Cpg ¼ 521:5 J/kgK Melting temperature Tm ¼ 1683:0 K Furnace pressure P0 ¼ 1333 Pa (10 Torr) Argon gas flow rate 0.025–1.0 Normal liters/min Chamber wall temperature 350.0 K Crucible inner diameter 0.072 m Crystal diameter 0.035 m Rotation rates of crystal and crucible 0 rpm Contact angle between melt and crystal yc ¼ 111 Emissivity e ¼ 0:55 (crystal, chamber wall), e ¼ 0:318 (melt), e ¼ 0:50 (quartz) e ¼ 0:9 (graphite, heater, pedestal, radiation shield, gas guide and puller)
3. Results A series of simulations under different gas flow rates has been conducted with and without the Marangoni effect and the gas guide. The effect of the gas guide with a simple shape will be discussed first. Figs. 2a and b show the isotherms and the stream function in gas phase with and without the
gas guide at 10 Torr and Q ¼ 0:5 l/min (l: 103 m3 at 0.1 MPa and 273 K). The gas guide changes the temperature and flow fields in the hot zone. The gas flow patterns in the upper chamber look similar. In both cases, a large toroidal roll cell caused by the buoyancy force occupies the empty chamber, and fresh gas flows down along the cold chamber wall. The flow pattern in the upper
3.20 3.61 3.20 3.61 Tmax;f ; Tmax;m : maximum temperatures in Cz furnaces and in the melt. DT: temperature difference between crucible wall and crystal wall across melt/gas interface. ½O max ; ½O min : maximum and minimum oxygen concentrations in the melt. ½O lc : oxygen concentration at the center of the melt/crystal interface. Vmax;m ; Vmax;s : maximum velocities in the melt and at the melt/gas interface. a
1.45 1016 4.86 1015 B B+Ma With a gas guide
16.36 16.30
1.003 1.054
1773.2 1771.5
1697.9 1696.8
14.7 13.8
2.52 1018 2.50 1018
1.54 1018 6.37 1016
2.17 5.03 5.04 5.05 4.39 1017 5.21 1017 5.21 1017 5.21 1017 2.57 1018 2.50 1018 2.50 1018 2.50 1018 16.4 13.1 13.1 13.2 1700.3 1696.8 1696.8 1696.9 1794.4 1790.9 1791.0 1791.1 1.008 1.425 1.012 0.573 17.27 17.15 17.16 17.16 B B+Ma Without gas guide
Heater power (kW) Driving mechanism Case
The contours of the stream function, the temperature and the oxygen concentration distribution in the melt, and the contours of SiO and the gas flow vectors above the hot zone are shown in Figs. 3 and 4. In Fig. 3, melt flow is driven only by the buoyancy force. Hereafter we denote this as the B-driven flow. Fig. 4 shows the result obtained with the Marangoni effect taken into account. We
Table 2 Summary of simulation resultsa
3.1. Hot zone in a conventional Cz furnace
Vs (mm/min)
gas chamber varies with the aspect ratio of the chamber. The present chamber has a larger aspect ratio than that is used by Bornside et al. [11]. They obtained a flow pattern characterized by the presence of two toroidal roll cells in the gas chamber. We obtained a similar flow pattern when we used the same aspect ratio as in Ref. [11]. The transport phenomena in these hot zones will be discussed in the following sections. The typical values in the simulations are listed in Table 2.
6.91 1017 5.79 1017 5.79 1017 5.79 1017
Vmax;m (cm/s) ½O lc (atoms/cm3) ½O min (atoms/cm3) ½O max (atoms/cm3) DT (K) Tmax;f (K)
Tmax;m (K)
Fig. 2. Isotherms over the global furnace (right) and the stream function in the gas phase (left). DT ¼ 72 K, DC ¼ 0:0005 m3/s, Cmax ¼ 0:0094 m3/s, Cmin ¼ 0:00053 m3/s. (a) Without gas guide. (b) With gas guide.
0.86 5.03 5.04 5.05
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46 Vmax;s (cm/s)
38
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
39
Fig. 3. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B-driven flow without a gas guide at a gas flow rate of Q ¼ 0:5 l/min, DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:0208 106 m3/s, Cmin ¼ 2:8487 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:8229; DySiO ¼ 0:02:
Fig. 4. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B+Ma-driven flow without a gas guide at a gas flow rate of Q ¼ 0:5 l/min. DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:0266 106 m3/s, Cmin ¼ 2:9547 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:7834; DySiO ¼ 0:02:
40
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
denote this situation as the B+Ma-driven flow. In both cases, there are two toroidal roll cells in the melt pool, one large cell flowing upward along the crucible wall and extending throughout most of the melt pool, and one small roll near the top of the crucible wall. It must be noted also that in these cases, the gas above the melt up to the limb of the crucible is almost stagnant. Fig. 3 indicates that the maximum surface velocity Vmax;s is low (see Table 2). This slow surface motion corresponds to the very low oxygen concentration on the surface. However, this lowoxygen-concentration melt is not carried to the crystal. Instead, the melt with much higher oxygen concentration ascends near the periphery of the crystal and transports oxygen to the melt/crystal interface. On the other hand, if the Marangoni effect is taken into account, radial surface flow is accelerated. Although the minimum oxygen concentration is slightly higher than that in Fig. 3, the low-oxygen-concentration melt on the surface is effectively carried to the area below the melt/ crystal interface. Then, the oxygen concentration in the crystal becomes lower than that in Fig. 3 (see Fig. 5). The accelerated melt flow increases the efficiency of heat transport from the crucible to the crystal. This enhanced heat transfer requires a reduction of the heater power, e.g., from 17.27 kW to 17.16 kW, as shown in Table 2. Simulation results with different heater powers, shown in Table 2, indicate that the pull rate is very sensitive to a very small change in heater power. However, it should be noted that the oxygen concentration at the center of crystal is quite insensitive to the pull rate. Fig. 5 shows the radial distribution of oxygen concentration (½O l atom/cm3) on the melt surface and at the melt/crystal interface, the interphase mass flux of oxygen or SiO (J [mol/m2 s]), and the mass-transfer coefficient kg defined as J ¼ kg cgi yiSiO ;
ð17Þ
where cgi is the molar density of gas at the surface temperature and yiSiO the molar fraction of SiO on the surface (yiSiO ¼ piSiO =P0 : piSiO is the partial pressure of SiO in equilibrium with xO (or ½O l ) at the surface temperature). These figures show that the oxygen concentration decreases very rapidly
Fig. 5. Radial distributions of oxygen concentration on the melt/crystal and at melt/gas interfaces, mass flux and gas-phase mass-transfer coefficient for cases without a gas guide at Q ¼ 0:5 l/min. (1) Radial distributions of [O]l (atoms/cm3) at the interfaces. (2) Radial distributions of mass flux J (mol/m2 s) of oxygen or SiO at the melt/gas interface. (3) Radial distributions of gas-phase mass-transfer coefficient, kg (m/s). Solid line: for Ra-driven flow; dotted lines: for B+Ma-driven flow.
near the crucible wall and remains almost constant over a wide area of the melt surface. However, neither the evaporation rate nor the evaporation coefficient is uniform along the surface. 3.2. Hot zone with a gas guide At the same gas flow rate (Q ¼ 0:5 l/min), the heater power is decreased by installing a cylindrical heat shield (gas guide). At the same time, the maximum temperature at the melt/crucible interface as well as the maximum oxygen concentration are both decreased, as shown in Table 2. Besides these thermal effects, the gas guide exerts significant effects on the melt flow pattern, as well as on the oxygen mass transfer and oxygen concentration, as shown in Figs. 6 and 7. The gas guide directs all of the gas to the melt surface. By flowing
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
41
Fig. 6. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B-driven flow with a gas guide at a gas flow rate of Q ¼ 0:5 l/min. DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:5094 106 m3/s, Cmin ¼ 2:1853 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:8016; DySiO ¼ 0:02:
Fig. 7. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B+Ma-driven flow with a gas guide at a gas flow rate of Q ¼ 0:5 l/min. DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:1352 106 m3/s, Cmin ¼ 2:2485 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:7908; DySiO ¼ 0:02:
42
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
along the melt surface, gas induces shear stress in the melt. In the B-driven case, this shear stress tends to initiate a secondary flow that is characterized by its outward radial surface velocity (see Fig. 6 and also Fig. 8 for a lower gas flow rate (Q ¼ 0:1 l/min)). Then, along the melt surface, the gas and melt concurrently flow over a wide area extending outward close to the crucible wall. This circulating flow creates a large pool of melt with a low oxygen concentration. With increasing gas flow rate, the area of the low-oxygen-concentration melt pool increases and the oxygen concentration decreases significantly (compare Fig. 6 with Fig. 8). However, the low-oxygen-concentration melt remains in the pool and is not efficiently conveyed to the melt/crystal interface. Highoxygen-concentration melt flows below the sheardriven circulation and ascends near the crystal periphery, carrying oxygen to the crystal/melt interface. Then, despite the low-oxygen-concentration melt on the surface, oxygen concentration remains high at the melt/crystal interface. In the B+Ma-driven flow case, there also appears a clockwise roll cell with low oxygen concentration in a restricted area below the gas guide, where the shear stress from gas flow becomes largest (see Fig. 7 and also Fig. 9 under Q ¼ 0:1 l/min). Due to the competition between the outward gas shear stress and the Marangoni effect acting toward the crystal on the surface, the flow pattern becomes more complex and the surface velocity becomes lower than that of the B+Ma-driven flow without a gas guide at the same gas flow rate (cf. Table 2). The surface is divided into three flow zones, i.e., (1) inward Ma-driven radial flow near the crystal, (2) outward shear-driven flow below the gas guide as a part of the shear-driven circulating flow, and (3) inward Ma-driven flow near the crucible wall. The low-oxygen-concentration melt is collected in the shear-driven circulating flow zone and the lowoxygen-concentration melt is efficiently conveyed to the area beneath the melt/crystal interface by the Marangoni-driven inward surface flow through a very thin layer near the crystal. Fig. 10 shows the radial distributions of oxygen concentration at the interface, interphase mass flux J and kg for B-driven and B+Ma-driven cases with a gas guide at a high gas flow rate, Q ¼ 0:5 l/
min. On comparing Fig. 10 with Fig. 5, it is clear that the gas guide enhances the gas phase mass transfer coefficient. In the B-driven case, negative mass flux, i.e., absorption of SiO, occurs at an area slightly away from the crucible wall. In this region, the thin flow of very low-oxygen-concentration melt submerges below the surface and the SiO concentration in the gas phase remains high because of the radial diffusion in the gas phase from the crucible wall. Then, part of the SiO once evaporated at the crucible wall is reabsorbed into the melt. A similar negative mass flux occurs in the B+Ma-driven case at the downstream end of the shear-driven circulating pool. At the boundaries of surface flow zones, kg shows significant change. 3.3. Effects of gas flow rate Figs. 11a and b show the effects of gas flow rate on the radial distribution of ½O l at melt/gas and melt/crystal interfaces over a wide range of gas flow rates from 0.025 to 0.75 l/min. On increasing the gas flow rate, oxygen concentration on the melt surface is reduced significantly in both B-driven and B+Ma-driven cases. However, oxygen concentration under the crystal is not the same as ½O l on the melt surface near the crystal and is strongly dependent on the melt flow pattern, as explained in the previous sections. Fig. 12 summarizes the effect of gas flow rate on ½O lc (the oxygen concentration at the melt/crystal interface on the axis (r ¼ 0)). At very low gas flow rates, ½O lc approaches the solubility limit, since SiO evaporation is highly suppressed due to the increased SiO concentration in the bulk of gas caused by poor ventilation. With increasing gas flow rate, the bulk concentration decreases and the driving force in the gas phase increases. Then, ½O lc decreases with increasing gas flow rate. At the same time, however, the increase of the gas flow rate is always accompanied by the increase of shear stress acting outward on the melt surface. In the B-driven case, the shear stress from the gas flow easily induces a clockwise circulating flow cell near the melt surface. Then, ½O lc remains almost constant after exhibiting a small minimum at around Q ¼ 0:05 l/min. On the other hand, in the B+Ma-driven case, the shear stress first reduces
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
43
Fig. 8. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B-driven flow with a gas guide at a gas flow rate of Q ¼ 0:1 l/min, DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:0492 106 m3/s, Cmin ¼ 2:2042 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:8148; DySiO ¼ 0:02:
Fig. 9. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B+Ma-driven flow with a gas guide at a gas flow rate of Q ¼ 0:1 l/min, DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:0180 106 m3/s, Cmin ¼ 2:3548 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:7704; DySiO ¼ 0:02:
[O]l (atoms/cm3)
44
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
1019 B
1018 1017
B+Ma
1016 1015
0
1
2
3
J (mol/m2s)
0.004 0.002 B 0.000
B+Ma
0
kg (m/s)
0.4
1
2
3 B
0.2 B+Ma
0.0 0
1
2
3
r (cm)
the inward melt surface flow. As the inward surface flow is suppressed, the concentrationboundary layer in the melt becomes thicker. As a result of reduced melt phase mass-transfer efficiency, the surface concentration decreases. At higher gas flow rates, gas shear initiates the sheardriven circulating flow under the gas guide and produces a rather small melt pool of low oxygen concentration. At any gas flow rate, the surface flows near the crucible and the crystal are governed by the Marangoni effect and always carry the surface melt to the area beneath the crystal. Thus, oxygen concentration on the melt surface, as well as ½O lc ; decreases as the gas flow rate increases if a weak Marangoni effect acts on the surface. However, further increase of gas flow rate beyond 0.5 l/min suppresses the Marangoni flow near the crystal and increases ½O lc again, as shown in Fig. 12. It must be noted that the
Fig. 11. Effect of gas flow rates on radial distributions of oxygen concentration at interfaces. Solid keys correspond to the results obtained without a gas guide. (a) Results for B-driven flow. (b) Results for B+Ma-driven flow.
1018 [O]lc (atoms/cm3)
Fig. 10. Radial distributions of oxygen concentration at the melt/crystal or melt/gas interfaces, mass flux and gas-phase mass-transfer coefficient with a gas guide at Q ¼ 0:5 l/min. (1) Radial distributions of [O]l (atoms/cm3) at the interfaces. (2) Radial distributions of mass flux J (mol/m2 s) of oxygen or SiO at the melt/gas interface. (3) Radial distributions of gas-phase mass-transfer coefficient, kg (m/s). Solid line: for B-driven flow. Dotted lines: for B+Ma-driven flow.
1017
10
16
0.0
B+Ma B B+Ma+Gas Guide B+Gas Guide 0.2
0.4
0.6
0.8
1.0
Q (L/min) Fig. 12. Effect of gas flow rate on the oxygen concentration at the center of the melt/crystal interface with a gas guide. Solid keys in the figure represent the results obtained without a gas guide.
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
magnitude of the Marangoni effect (i.e., temperature coefficient of surface tension, gT ) is an important parameter to determine the melt flow pattern and ultimately the oxygen concentration in this small Cz furnace.
4. Discussion The present global analyses reveal significant effects of gas-phase transport phenomena and the Marangoni effect on thermal environment, flow pattern, oxygen transport and oxygen concentration at the melt/crystal interface. However, these analyses have been conducted for a small silicon Cz furnace without rotation of the crystal and the crucible. The results may not be valid in larger Cz furnaces. The simulations will be further extended in order to include crystal and crucible rotations. Further development of the numerical code is required to enable simulations of much larger industrial Cz furnaces. In order to realize global analysis in large furnaces, we must develop a more efficient numerical code for calculating threedimensional time-dependent flow in the melt and turbulent flow in the gas phase, as well as a robust algorithm for determining time-dependent threedimensional melt/crystal interface shape.
45
flow and oxygen transport. If there were no Marangoni effect, the increase of the gas flow rate would increase the oxygen concentration in the crystal. If the Marangoni effect worked on the surface, oxygen concentration would decrease significantly by increasing the gas flow rate. These results indicate that gas flow must always be taken into account in the analysis of a Cz furnace with a heat shield. 3. Global analyses revealed that the evaporation rate of SiO is not uniform over the melt surface. In some cases, the reabsorption of SiO occurs near the crucible wall. Further simulation and development of numerical codes are needed to extend the global analyses to more realistic conditions, such as those with crucible and crystal rotations in large industrial Cz furnaces. Acknowledgements This work was supported by JSPS Research for the Future Program in the field of Atomic Scale Surface and Interface Dynamics. M.W. Li expresses his gratitude for JSPS’s support as a Postdoctoral Researcher at Kyushu University. References
5. Conclusions A global simulation code was developed to simulate the complex transport phenomena in a small silicon Cz furnace. The code could treat pseudosteady axisymmetric systems. The results revealed the effects of gas flow and the Marangoni effect on the melt flow pattern and the mass transfer characteristics of oxygen impurity. 1. The Marangoni effect (evaluated using a small temperature coefficient of surface tension, gT ¼ 7 105 N/mK) induced marked effects on the melt flow and oxygen concentration in a small Cz furnace. 2. A cylindrical heat shield, which acted as a gas guide, caused a significant decrease of heater power. It also exerted significant effect on gas
[1] N. Kobayashi, J. Crystal Growth 108 (1991) 240. [2] P. Sabhapathy, M.E. Salcudean, J. Crystal Growth 113 (1991) 164. [3] H. Hirata, K. Hoshikawa, J. Crystal Growth 125 (1992) 181. [4] T.A. Kinney, R.A. Brown, J. Crystal Growth 132 (1993) 551. [5] K. Kakimoto, K.W. Yi, M. Eguchi, J. Crystal Growth 163 (1996) 238. [6] K.W. Yi, K. Kakimoto, M. Eguchi, H. Noguchi, J. Crystal Growth 165 (1996) 358. [7] M. Watanabe, K.W. Yi, T. Hibiya, K. Kakimoto, Prog. Crystal Growth Characterization Mater. 38 (1999) 215. [8] S. Kawanishi, S. Togawa, K. Izunome, K. Terashima, S. Kimura, Jpn. J. Appl. Phys. 34 (1995) 5885. [9] I. Kanda, T. Suzuki, K. Kojima, J. Crystal Growth 166 (1996) 669. [10] N. Machida, Y. Suzuki, K. Abe, N. Ono, M. Kida, Y. Shimizu, J. Crystal Growth 186 (1998) 362. [11] D.E. Bornside, R.A. Brown, T. Fujiwara, H. Fujiwara, T. Kubo, J. Electrochem. Soc. 142 (1995) 2790.
46
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
[12] T. Tsukada, M. Hozawa, N. Imaishi, J. Chem. Eng. Japan 27 (1994) 25. [13] H. Honda, N. Imaishi, T. Tsukada, M. Hozawa, J. Chem. Eng. Japan 25 (1992) 84. [14] A. Muhe, R. Backofen, J. Fainberg, G. Muller, E. Dornberger, E. Tomzig, W. von Ammon, J. Crystal Growth 198/199 (1999) 409.
[15] T. Carlberg, J. Electrochem. Soc. 133 (1986) 1940. [16] M. Ratto, E. Ricci, E. Arato, J. Crystal Growth 217 (2000) 233. [17] K. Nogi, private communication. [18] H. Fujii, T. Matsumoto, N. Hata, M. Kohno, K. Nogi, Metall. Mater. Trans. A 31 (2000) 1585.
Global simulation of a silicon Czochralski furnace Mingwei Lia, Yourong Lia, Nobuyuki Imaishia,*, Takao Tsukadab a
Institute of Advanced Material Study, Kyushu University, 6-1 Kasuga-koen, Kasuga, Fukuoka 816-8580, Japan Institute of Chemical Reaction Science, Tohoku University, 2-1-1 Katahira, Aoba-ku, Sendai 980-8577, Japan
b
Received 6 July 2000; accepted 25 July 2001 Communicated by J.J. Derby
Abstract To understand the characteristics of the Czochralski (Cz) furnace for the single-crystal growth of silicon, a set of global analyses of momentum, heat and mass transfer in small Cz furnaces (crucible diameter: 7.2 cm, crystal diameter: 3.5 cm, operated in a 10 Torr argon flow environment) is carried out using the finite-element method. The global analysis assumes a pseudosteady axisymmetric state with laminar flow, equilibrium relations at the melt/silica interface and vapor–liquid chemical equilibrium at the melt/gas interface, as well as the segregation coefficient of unity at the melt/crystal interface. Convective and conductive heat transfers, radiative heat transfer between diffuse surfaces and the Navier–Stokes equations for gas and melt phases are all combined and solved together. Thus, the velocities and temperatures obtained are used to calculate the oxygen concentrations. The global analysis code is effectively used to discuss the influences of the Marangoni effect and a gas guide (or a heat shield) placed between the crucible and the crystal. The results indicate that the gas guide reduces the heater power and changes the melt flow pattern and oxygen transport. The melt flow pattern is strongly dependent on the Marangoni effect and gas flow near the surface, and changes the oxygen concentration significantly. This analysis reveals the importance and effectiveness of global analysis. r 2002 Elsevier Science B.V. All rights reserved. Keywords: A1. Computer simulation; A1. Fluid flows; A1. Heat transfer; A1. Mass transfer; A2. Czochralski method; B2. Semiconducting silicon
1. Introduction For the production of high-quality silicon crystals by the Czochralski (Cz) method, the control of concentration and uniformity of oxygen, as well as of the temperature gradient, in the growing crystal is required. In order to control these factors, most of the recent Cz furnaces are equipped with substructures in the furnace, such as *Corresponding author. Tel: +81-92-583-7793; fax: +81-92583-7796. E-mail address: [email protected] (N. Imaishi).
a funnel-shaped radiation shield and multiple heaters for the sake of thermal environment control as well as for saving energy. To control oxygen concentration, the operating pressure and the gas flow rate are included as operating parameters, in addition to conventional control factors such as crystal and crucible rotation rates. The application of various types of magnetic fields to the electrically conducting silicon melt is fast becoming popular for stabilizing turbulent melt flow and controlling oxygen transfer in large-scale crucibles. The oxygen in the crystal originates from the dissolution of silica from the crucible into
0022-0248/02/$ - see front matter r 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 0 1 ) 0 1 6 3 4 - 7
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
molten silicon. Dissolved oxygen is convectively transported to the melt/crystal interface and incorporated into the crystal. It is believed that before arriving at the melt/crystal interface, most of the dissolved oxygen is evaporated into the gas phase through the melt/gas interface via the chemical reaction SiðlÞ þ OðlÞ"SiOðgÞ at the melt/gas interface. In recent years, many numerical analyses [1–7] have been conducted on melt flow and oxygen transport in the silicon melt. These analyses have helped us greatly in understanding the transport mechanism. The experimental results obtained by Kawanishi et al. [8] and Kanda et al. [9] suggest that the melt surface flow plays a very important role in the radial nonuniformity of oxygen concentration in crystal. However, the boundary condition for the diffusion equation at the melt/gas interface in these numerical simulations assumed either a zero oxygen concentration [1,3,4] or a constant evaporation coefficient [5–7]. The zero oxygen concentration boundary condition implies that there is no diffusion resistance in the gas phase. This is in contrast to the experimental results obtained by Machida et al. [10], which revealed that the gas phase pressure and the gas flow rate significantly affect the oxygen concentration in the crystal. In their experiments, the increase of gas pressure or decrease of gas flow rate caused a considerable increase of oxygen concentration in the crystal. In a crystal grown in a Cz furnace with a heat shield, which also acts as a gas guide, the oxygen concentration decreases until a critical gas flow rate is reached. Beyond the critical flow rate, however, the oxygen concentration increases again with further increase of the gas flow rate. Machida et al. speculated that at higher gas flow rates, the shear stress from the gas flow may change the melt flow pattern, decrease the evaporation rate of SiO and consequently increase the oxygen concentration in the grown crystal. Bornside et al. [11] conducted the first numerical simulation of the evaporation of SiO and the incorporation of carbon through the gas phase in a conventional Cz furnace. However, in their analysis, they assumed that the melt surface was rigid with a constant vapor pressure of SiO. The effect of substructures in the Cz furnace, such as a heat shield, has not yet been analyzed,
33
except for some analyses on their effect on the melt/crystal interface shape [12,13]. The basic structure of the Cz furnace is very simple and the number of components is rather small. Despite such a simple structure, phenomena in the Cz furnace are highly nonlinear and they are strongly coupled with each other. Therefore, one small change may cause various effects. For example, a cylindrical or funnel-shaped heat shield changes the temperature gradient in the crystal as well as the shape of the melt/crystal interface [12,13]. At the same time, the heat shield changes the path of gas flow near the melt/gas interface as well as the evaporation rate of oxygen. The heat shield reduces the radiative heat loss and in turn, the heater power, temperature distributions and melt flow are also modified. These change the solubility of SiO2 in silicon melt, the equilibrium vapor pressure of SiO at the melt/gas interface, as well as the gas-phase mass-transfer coefficient. These are all closely related to each other and no individual component analysis can claim its own significance. In order to understand the overall characteristics of the Cz furnace, ‘‘global analysis’’ must be applied. Previously, Kinney and Brown [4] and Muhe et al. [14] analyzed the melt convection, and conductive and radiative heat transfer in a Cz furnace, but did not pay any attention to the gasphase transport phenomena. In the present work, a series of simulations was conducted for a small silicon Cz furnace using a laboratory-made global analysis code that includes heat exchange in the entire furnace, melt and gas flows, and oxygen and SiO transport in the melt and in the gas. The code requires input data for only one set of operating conditions, such as gas flow rate, gas-phase pressure, heater power, furnace chamber wall temperature, geometry of the furnace, diameter and length of the crystal, and melt volume. The pseudosteady state and axisymmetry were assumed. The melt/crystal interface shape was determined as the isotherm at the melting temperature that satisfies the heat balance as well as all boundary conditions. The effects of a gas guide with a simple shape were discussed, together with the Marangoni (thermocapillary) effect. For simplicity, crystal and crucible rotations were neglected in this work.
34
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
2. Model formulation 2.1. Basic assumptions The following assumptions are introduced to the model: (1) The flow is laminar in both the melt and the gas and the melt is incompressible. (2) Thermophysical properties of the melt are assumed to be constant except for the temperature dependence of density and surface tension, in the calculation of the buoyancy force and the Marangoni effect, respectively. However, the thermophysical properties of argon gas are treated as temperature-dependent. (3) The PVT relation of argon gas is expressed by the ideal gas law. (4) No supercooling occurs at the melt/crystal interface and the interface shape coincides with the melting temperature isotherm. (5) Oxygen concentration in the argon gas at the inlet is zero. (6) Oxygen concentration at the melt/crucible interface is given by the solubility equilibrium [15]. (7) The oxygen concentration and the partial pressure of SiO are in chemical equilibrium at the melt–gas interface [11]. (8) The partial pressure of oxygen in the bulk of argon gas is zero. (9) The segregation coefficient of oxygen is assumed to be unity. (10) The temperature coefficient of surface tension is taken as a disposable parameter to test the role of the Marangoni effect. (11) Evaporated SiO does not deposit on any solid surface. Among these assumptions, only the last one is not realistic but may be acceptable as long as the discussion is limited to the evaporation rate of SiO. For further analysis, the vapor pressure of SiO and its reaction with carbon materials should be included, as in Ref. [11]. Assumption (9) is introduced here because there is no reliable chemical analysis data on the argon gas used in the experiments. If such data were available, the oxygen supply from the gas phase must be taken into account, as in Ref. [16]. 2.2. Mathematical models 2.2.1. Governing equations Under the above assumptions, the general governing equations for flow, temperature and concentration fields are given as follows.
In the melt: r vl ¼ 0;
ð1Þ
rl vl rvl ¼ rpl r sl þ rl gbðTl Tm Þez ;
ð2Þ
rl Cpl vl rTl ¼ ll r rTl ;
ð3Þ
vl rxO ¼ DO r rxO :
ð4Þ
In the gas: r ðrg vg Þ ¼ 0;
ð5Þ
rg vg rvg ¼ rpg r sg rg gez ;
ð6Þ
rg Cpg vg rTg ¼ r ðkg rTg Þ;
ð7Þ
r ðcg vg ySiO Þ ¼ r ðcg DSiO rySiO Þ:
ð8Þ
In the crystal: rs Cps Vs ez rTs ¼ r ðls rTs Þ:
ð9Þ
In the heater: r ðlh rTh Þ þ Q ¼ 0:
ð10Þ
In the other solid materials: r ðli rTi Þ ¼ 0:
ð11Þ
T; v; p; and s are the temperature, velocity vector, pressure and stress tensor, respectively. Tm is the melting temperature, r the density, l the thermal conductivity, D the mass diffusivity, b the thermal expansion coefficient, g the gravitational acceleration constant, Cp the heat capacity, Vs the crystal pulling rate and c the molar density. x and y are the mole fractions of oxygen and SiO, respectively. Subscripts l; g; s; h; i; O and SiO indicate melt, gas, crystal, heater, other materials, oxygen and SiO, respectively. ez is a unit vector parallel to the z axis. 2.2.2. Boundary conditions for the temperature field At the melt/gas interface: ll n rTl ¼ lg n rTg þ qrad;l :
ð12aÞ
At the melt/crystal interface: ll n rTl ls n rTs ¼ rs Vs DHs n ez ;
ð12bÞ
Tl ¼ Ts ¼ Tm :
ð12cÞ
35
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
At the crystal side and top: ls n rTs ¼ lg n rTg þ qrad;s : At the other exposed surfaces: li n rTi ¼ lg n rTg þ qrad;i :
ð12eÞ
At the inlet of Ar gas: Tg ¼ Tg0 :
ð12fÞ
At the outlet of gas: n rTg ¼ 0:
tl : ney tg : ney ¼ gT rTl ey ;
ð14kÞ
ey vl ey vg ¼ 0:
ð14lÞ
ð12dÞ
ð12gÞ
In these equations, qrad;i is the net radiative heat flux per unit area on the surfaces, and is given by " ! # X 4 4 qrad;i ¼ ei sB Ti Aj ej Tj Gji =Ai =ei ; ð13Þ
At the melt/crystal interface: u1 ¼ 0;
ð14mÞ
v1 ¼ 0;
ð14nÞ
wl ¼ rOs :
ð14oÞ
At the other surfaces of crystal and the puller surface: ug ¼ 0;
ð14pÞ
vg ¼ 0;
ð14qÞ
wg ¼ rOs :
ð14rÞ
j
where Gji is Gebhart’s absorption factor, which is the fraction of emission from surface Aj that reaches Ai and is absorbed. Gji as well as the view factors between any two surface elements in the furnace is calculated using our code [13]. DHs is the latent heat of fusion, n the unit normal vector of the interface element and e the emissivity. 2.2.3. Boundary conditions for the flow field At the side wall and the bottom of the crucible:
At the inlet of Ar gas: ug ¼ 0;
ð14sÞ
wg ¼ 0;
ð14tÞ
vg ¼ Vin :
ð14uÞ
At the outlet of Ar gas: t n ¼ 0:
ð14vÞ
u1 ¼ 0;
ð14aÞ
v1 ¼ 0;
ð14bÞ
ug ¼ 0;
ð14wÞ
wl ¼ rOc :
ð14cÞ
vg ¼ 0;
ð14xÞ
At the other crucible surfaces and the pedestal surface:
wg ¼ 0:
ð14yÞ
ug ¼ 0;
ð14dÞ
vg ¼ 0;
ð14eÞ
wg ¼ rOc :
ð14fÞ
u; v; and w are the radial, axial and the azimuthal components of velocity, and Os and Oc are the angular rotation rates of the crystal and the crucible, respectively. t and ey are the unit tangent and azimuthal vectors of each interface element, respectively. The melt/crystal interface shape is determined such that Eq. (12c) is satisfied, i.e., the interface coincides with the melting-point isotherm.
At the melt/gas interface:
At other motionless solid surfaces:
n vl ¼ 0;
ð14gÞ
n vg ¼ 0;
ð14hÞ
tl : nt tg : nt ¼ gT rTl t;
ð14iÞ
2.2.4. Boundary conditions for oxygen and SiO At the crucible wall [15]:
t vl t vg ¼ 0;
ð14jÞ
xO ¼ 0:01eð7:822 370=TÞ :
ð15aÞ
36
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
At the melt/gas interface [11]: ySiO ¼ 101 325xO e
ð17:821 000:0=TÞ
=P0 ;
cg DSiO n rySiO ¼ cl DO n rxO :
2.4. Simulation conditions ð15bÞ ð15cÞ
At the melt/crystal interface: n rxO ¼ 0:
ð15dÞ
At the inlet of Ar gas: ySiO ¼ 0:
ð15eÞ
At the outlet of gas: qySiO ¼ 0: qz At the other exposed surfaces:
ð15fÞ
n rySiO ¼ 0;
ð15gÞ
where P0 (=1333 Pa or 10 Torr) is the total pressure. 2.3. Numerical procedure To solve this problem, the Galerkin finiteelement method is used. The calculation domain is discretized by isoparametric quadrilateral elements. In each element, velocity vectors, temperature and concentration are approximated with the bilinear trial function f as follows: X vðr; zÞ ¼ f i vi ; ð16aÞ X ð16bÞ Tðr; zÞ ¼ fi Ti ; X ð16cÞ yðr; zÞ ¼ fi ySiOi ; X ð16dÞ xðr; zÞ ¼ fi xOi : The pressure is assumed to be constant in each element. By substituting Eqs. (16a)–(16d) into the governing equations and the boundary conditions for the flow, temperature and concentration fields, the Galerkin procedure gives a set of algebraic equations. The melt/gas interface shape is calculated by solving the Young–Laplace equation using the Galerkin method with the Hermit cubic trial function. In the following analyses, this melt/gas interface shape is approximated by 50 linear segments.
The configurations of the Cz furnaces adopted in the present global analyses are shown in Fig. 1. The inner diameter of both crucibles is 7.2 cm, and the diameter of both crystals is 3.5 cm. The heat shield (or gas guide) is a straight cylinder (inner diameter 5.15 cm, thickness 2.0 mm) and its lower end is located 1.14 cm above the melt surface, as shown in Fig. 1. It is unrealistic to install a heat shield in such a small furnace. The heat shield is adopted in this analysis to elucidate the distinguishable effects of such a substructure on heat transfer and mass transfer, and also to demonstrate the importance of the global analysis. In this paper, crystal and crucible rotations, as well as any azimuthal distribution, are neglected. Thus, we assume Os ¼ Oc ¼ 0 and neglect Eqs. (14k) and (14l). In the present work, the diameter and length of the crystal and the heater power are given. The crystal pull rate is treated as one of the unknown
Fig. 1. Configuration of the Cz furnaces used in this study.
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
variables, together with velocities, pressures, temperatures, and the melt/crystal interface coordinates, by simultaneously solving the set of nonlinear algebraic equations using the Newton– Raphson method. The total number of elements, nodal points and unknowns are 54 542, 55 107 and 242 486, respectively. The concentration field is computed by using much finer elements, especially in the melt phase, with 81 576 elements, 82 458 nodal points and 82 458 unknown variables. The sizes of the minimum elements that are located at the corner of the melt/gas interface and the crucible wall in the melt are Drmin ¼ 0:125 mm, Dzmin ¼ 0:073 mm for the flow and heat transfer analysis, and Drmin ¼ 2:60 mm, Dzmin ¼ 0:71 mm for the oxygen transport analysis. For mass transfer calculations, the temperature and velocities at each point are interpolated based on the obtained velocity and temperature values. The code requires only a set of operating conditions as input data, such as gas flow rate, gas phase pressure, heater power, furnace chamber wall temperature, geometry of the furnace, diameter and length of the crystal, and total volume of silicon (melt and crystal). Typical conditions are listed in Table 1, together with the thermophysical properties used in this work. The temperature coefficient of surface tension is assumed to be either zero or 7 105 N/mK. The latter value is taken from the recent experiments of Nogi [17] (In their final report [18], the value is revised to – 6.2 105 N/mK.). The heater power distribution is evaluated by solving the Laplace equation for the current in the heater, and then the power distribution is averaged over both azimuthal and radial directions expressed as a function of the longitudinal coordinate.
37
Table 1 Physical properties and processing parameters Thermal conductivity l ¼ 22:0 W/mK (crystal), l ¼ 64U0 W/mK (melt), l ¼ 2:98 W/mK (quartz) l ¼ 60:0 W/mK (graphite, heater, pedestal, radiation shield and gas guide) l ¼ 59:0 W/mK (puller) Viscosity m ¼ 7U0 104 kg/ms (melt) Diffusivity DO ¼ 5 108 m2/s, DSiO ¼ 8:626 106 T 1:75 =P0 m2/s Density r ¼ 2530:0 kg/m3 (melt), r ¼ P0 M=RT kg/m3 (Argon gas) Thermal expansion coefficient of melt bl ¼ 1:5 104 K1 Surface tension and its temperature coefficient g ¼ 0:735 N/m, gT ¼ 0:0 or 7.0 105 N/mK Heat of fusion DHm ¼ 1410 J/kg Heat capacity Cpl ¼ 1000 J/kgK, Cps ¼ 1000 J/kgK, Cpg ¼ 521:5 J/kgK Melting temperature Tm ¼ 1683:0 K Furnace pressure P0 ¼ 1333 Pa (10 Torr) Argon gas flow rate 0.025–1.0 Normal liters/min Chamber wall temperature 350.0 K Crucible inner diameter 0.072 m Crystal diameter 0.035 m Rotation rates of crystal and crucible 0 rpm Contact angle between melt and crystal yc ¼ 111 Emissivity e ¼ 0:55 (crystal, chamber wall), e ¼ 0:318 (melt), e ¼ 0:50 (quartz) e ¼ 0:9 (graphite, heater, pedestal, radiation shield, gas guide and puller)
3. Results A series of simulations under different gas flow rates has been conducted with and without the Marangoni effect and the gas guide. The effect of the gas guide with a simple shape will be discussed first. Figs. 2a and b show the isotherms and the stream function in gas phase with and without the
gas guide at 10 Torr and Q ¼ 0:5 l/min (l: 103 m3 at 0.1 MPa and 273 K). The gas guide changes the temperature and flow fields in the hot zone. The gas flow patterns in the upper chamber look similar. In both cases, a large toroidal roll cell caused by the buoyancy force occupies the empty chamber, and fresh gas flows down along the cold chamber wall. The flow pattern in the upper
3.20 3.61 3.20 3.61 Tmax;f ; Tmax;m : maximum temperatures in Cz furnaces and in the melt. DT: temperature difference between crucible wall and crystal wall across melt/gas interface. ½O max ; ½O min : maximum and minimum oxygen concentrations in the melt. ½O lc : oxygen concentration at the center of the melt/crystal interface. Vmax;m ; Vmax;s : maximum velocities in the melt and at the melt/gas interface. a
1.45 1016 4.86 1015 B B+Ma With a gas guide
16.36 16.30
1.003 1.054
1773.2 1771.5
1697.9 1696.8
14.7 13.8
2.52 1018 2.50 1018
1.54 1018 6.37 1016
2.17 5.03 5.04 5.05 4.39 1017 5.21 1017 5.21 1017 5.21 1017 2.57 1018 2.50 1018 2.50 1018 2.50 1018 16.4 13.1 13.1 13.2 1700.3 1696.8 1696.8 1696.9 1794.4 1790.9 1791.0 1791.1 1.008 1.425 1.012 0.573 17.27 17.15 17.16 17.16 B B+Ma Without gas guide
Heater power (kW) Driving mechanism Case
The contours of the stream function, the temperature and the oxygen concentration distribution in the melt, and the contours of SiO and the gas flow vectors above the hot zone are shown in Figs. 3 and 4. In Fig. 3, melt flow is driven only by the buoyancy force. Hereafter we denote this as the B-driven flow. Fig. 4 shows the result obtained with the Marangoni effect taken into account. We
Table 2 Summary of simulation resultsa
3.1. Hot zone in a conventional Cz furnace
Vs (mm/min)
gas chamber varies with the aspect ratio of the chamber. The present chamber has a larger aspect ratio than that is used by Bornside et al. [11]. They obtained a flow pattern characterized by the presence of two toroidal roll cells in the gas chamber. We obtained a similar flow pattern when we used the same aspect ratio as in Ref. [11]. The transport phenomena in these hot zones will be discussed in the following sections. The typical values in the simulations are listed in Table 2.
6.91 1017 5.79 1017 5.79 1017 5.79 1017
Vmax;m (cm/s) ½O lc (atoms/cm3) ½O min (atoms/cm3) ½O max (atoms/cm3) DT (K) Tmax;f (K)
Tmax;m (K)
Fig. 2. Isotherms over the global furnace (right) and the stream function in the gas phase (left). DT ¼ 72 K, DC ¼ 0:0005 m3/s, Cmax ¼ 0:0094 m3/s, Cmin ¼ 0:00053 m3/s. (a) Without gas guide. (b) With gas guide.
0.86 5.03 5.04 5.05
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46 Vmax;s (cm/s)
38
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
39
Fig. 3. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B-driven flow without a gas guide at a gas flow rate of Q ¼ 0:5 l/min, DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:0208 106 m3/s, Cmin ¼ 2:8487 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:8229; DySiO ¼ 0:02:
Fig. 4. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B+Ma-driven flow without a gas guide at a gas flow rate of Q ¼ 0:5 l/min. DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:0266 106 m3/s, Cmin ¼ 2:9547 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:7834; DySiO ¼ 0:02:
40
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
denote this situation as the B+Ma-driven flow. In both cases, there are two toroidal roll cells in the melt pool, one large cell flowing upward along the crucible wall and extending throughout most of the melt pool, and one small roll near the top of the crucible wall. It must be noted also that in these cases, the gas above the melt up to the limb of the crucible is almost stagnant. Fig. 3 indicates that the maximum surface velocity Vmax;s is low (see Table 2). This slow surface motion corresponds to the very low oxygen concentration on the surface. However, this lowoxygen-concentration melt is not carried to the crystal. Instead, the melt with much higher oxygen concentration ascends near the periphery of the crystal and transports oxygen to the melt/crystal interface. On the other hand, if the Marangoni effect is taken into account, radial surface flow is accelerated. Although the minimum oxygen concentration is slightly higher than that in Fig. 3, the low-oxygen-concentration melt on the surface is effectively carried to the area below the melt/ crystal interface. Then, the oxygen concentration in the crystal becomes lower than that in Fig. 3 (see Fig. 5). The accelerated melt flow increases the efficiency of heat transport from the crucible to the crystal. This enhanced heat transfer requires a reduction of the heater power, e.g., from 17.27 kW to 17.16 kW, as shown in Table 2. Simulation results with different heater powers, shown in Table 2, indicate that the pull rate is very sensitive to a very small change in heater power. However, it should be noted that the oxygen concentration at the center of crystal is quite insensitive to the pull rate. Fig. 5 shows the radial distribution of oxygen concentration (½O l atom/cm3) on the melt surface and at the melt/crystal interface, the interphase mass flux of oxygen or SiO (J [mol/m2 s]), and the mass-transfer coefficient kg defined as J ¼ kg cgi yiSiO ;
ð17Þ
where cgi is the molar density of gas at the surface temperature and yiSiO the molar fraction of SiO on the surface (yiSiO ¼ piSiO =P0 : piSiO is the partial pressure of SiO in equilibrium with xO (or ½O l ) at the surface temperature). These figures show that the oxygen concentration decreases very rapidly
Fig. 5. Radial distributions of oxygen concentration on the melt/crystal and at melt/gas interfaces, mass flux and gas-phase mass-transfer coefficient for cases without a gas guide at Q ¼ 0:5 l/min. (1) Radial distributions of [O]l (atoms/cm3) at the interfaces. (2) Radial distributions of mass flux J (mol/m2 s) of oxygen or SiO at the melt/gas interface. (3) Radial distributions of gas-phase mass-transfer coefficient, kg (m/s). Solid line: for Ra-driven flow; dotted lines: for B+Ma-driven flow.
near the crucible wall and remains almost constant over a wide area of the melt surface. However, neither the evaporation rate nor the evaporation coefficient is uniform along the surface. 3.2. Hot zone with a gas guide At the same gas flow rate (Q ¼ 0:5 l/min), the heater power is decreased by installing a cylindrical heat shield (gas guide). At the same time, the maximum temperature at the melt/crucible interface as well as the maximum oxygen concentration are both decreased, as shown in Table 2. Besides these thermal effects, the gas guide exerts significant effects on the melt flow pattern, as well as on the oxygen mass transfer and oxygen concentration, as shown in Figs. 6 and 7. The gas guide directs all of the gas to the melt surface. By flowing
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
41
Fig. 6. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B-driven flow with a gas guide at a gas flow rate of Q ¼ 0:5 l/min. DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:5094 106 m3/s, Cmin ¼ 2:1853 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:8016; DySiO ¼ 0:02:
Fig. 7. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B+Ma-driven flow with a gas guide at a gas flow rate of Q ¼ 0:5 l/min. DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:1352 106 m3/s, Cmin ¼ 2:2485 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:7908; DySiO ¼ 0:02:
42
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
along the melt surface, gas induces shear stress in the melt. In the B-driven case, this shear stress tends to initiate a secondary flow that is characterized by its outward radial surface velocity (see Fig. 6 and also Fig. 8 for a lower gas flow rate (Q ¼ 0:1 l/min)). Then, along the melt surface, the gas and melt concurrently flow over a wide area extending outward close to the crucible wall. This circulating flow creates a large pool of melt with a low oxygen concentration. With increasing gas flow rate, the area of the low-oxygen-concentration melt pool increases and the oxygen concentration decreases significantly (compare Fig. 6 with Fig. 8). However, the low-oxygen-concentration melt remains in the pool and is not efficiently conveyed to the melt/crystal interface. Highoxygen-concentration melt flows below the sheardriven circulation and ascends near the crystal periphery, carrying oxygen to the crystal/melt interface. Then, despite the low-oxygen-concentration melt on the surface, oxygen concentration remains high at the melt/crystal interface. In the B+Ma-driven flow case, there also appears a clockwise roll cell with low oxygen concentration in a restricted area below the gas guide, where the shear stress from gas flow becomes largest (see Fig. 7 and also Fig. 9 under Q ¼ 0:1 l/min). Due to the competition between the outward gas shear stress and the Marangoni effect acting toward the crystal on the surface, the flow pattern becomes more complex and the surface velocity becomes lower than that of the B+Ma-driven flow without a gas guide at the same gas flow rate (cf. Table 2). The surface is divided into three flow zones, i.e., (1) inward Ma-driven radial flow near the crystal, (2) outward shear-driven flow below the gas guide as a part of the shear-driven circulating flow, and (3) inward Ma-driven flow near the crucible wall. The low-oxygen-concentration melt is collected in the shear-driven circulating flow zone and the lowoxygen-concentration melt is efficiently conveyed to the area beneath the melt/crystal interface by the Marangoni-driven inward surface flow through a very thin layer near the crystal. Fig. 10 shows the radial distributions of oxygen concentration at the interface, interphase mass flux J and kg for B-driven and B+Ma-driven cases with a gas guide at a high gas flow rate, Q ¼ 0:5 l/
min. On comparing Fig. 10 with Fig. 5, it is clear that the gas guide enhances the gas phase mass transfer coefficient. In the B-driven case, negative mass flux, i.e., absorption of SiO, occurs at an area slightly away from the crucible wall. In this region, the thin flow of very low-oxygen-concentration melt submerges below the surface and the SiO concentration in the gas phase remains high because of the radial diffusion in the gas phase from the crucible wall. Then, part of the SiO once evaporated at the crucible wall is reabsorbed into the melt. A similar negative mass flux occurs in the B+Ma-driven case at the downstream end of the shear-driven circulating pool. At the boundaries of surface flow zones, kg shows significant change. 3.3. Effects of gas flow rate Figs. 11a and b show the effects of gas flow rate on the radial distribution of ½O l at melt/gas and melt/crystal interfaces over a wide range of gas flow rates from 0.025 to 0.75 l/min. On increasing the gas flow rate, oxygen concentration on the melt surface is reduced significantly in both B-driven and B+Ma-driven cases. However, oxygen concentration under the crystal is not the same as ½O l on the melt surface near the crystal and is strongly dependent on the melt flow pattern, as explained in the previous sections. Fig. 12 summarizes the effect of gas flow rate on ½O lc (the oxygen concentration at the melt/crystal interface on the axis (r ¼ 0)). At very low gas flow rates, ½O lc approaches the solubility limit, since SiO evaporation is highly suppressed due to the increased SiO concentration in the bulk of gas caused by poor ventilation. With increasing gas flow rate, the bulk concentration decreases and the driving force in the gas phase increases. Then, ½O lc decreases with increasing gas flow rate. At the same time, however, the increase of the gas flow rate is always accompanied by the increase of shear stress acting outward on the melt surface. In the B-driven case, the shear stress from the gas flow easily induces a clockwise circulating flow cell near the melt surface. Then, ½O lc remains almost constant after exhibiting a small minimum at around Q ¼ 0:05 l/min. On the other hand, in the B+Ma-driven case, the shear stress first reduces
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
43
Fig. 8. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B-driven flow with a gas guide at a gas flow rate of Q ¼ 0:1 l/min, DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:0492 106 m3/s, Cmin ¼ 2:2042 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:8148; DySiO ¼ 0:02:
Fig. 9. Contours of stream function, C; temperature (color scale on the right) and oxygen concentration (xO ; color scale on the left) in the melt, and the contours of ySiO and the gas flow vectors above the melt/gas interface in the case of B+Ma-driven flow with a gas guide at a gas flow rate of Q ¼ 0:1 l/min, DxO ¼ 2:56 1017 atoms/cm3, Cmax ¼ 0:0180 106 m3/s, Cmin ¼ 2:3548 106 m3/s, DC ¼ 0:07873 106 m3/s, ySiO;max ¼ 0:7704; DySiO ¼ 0:02:
[O]l (atoms/cm3)
44
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
1019 B
1018 1017
B+Ma
1016 1015
0
1
2
3
J (mol/m2s)
0.004 0.002 B 0.000
B+Ma
0
kg (m/s)
0.4
1
2
3 B
0.2 B+Ma
0.0 0
1
2
3
r (cm)
the inward melt surface flow. As the inward surface flow is suppressed, the concentrationboundary layer in the melt becomes thicker. As a result of reduced melt phase mass-transfer efficiency, the surface concentration decreases. At higher gas flow rates, gas shear initiates the sheardriven circulating flow under the gas guide and produces a rather small melt pool of low oxygen concentration. At any gas flow rate, the surface flows near the crucible and the crystal are governed by the Marangoni effect and always carry the surface melt to the area beneath the crystal. Thus, oxygen concentration on the melt surface, as well as ½O lc ; decreases as the gas flow rate increases if a weak Marangoni effect acts on the surface. However, further increase of gas flow rate beyond 0.5 l/min suppresses the Marangoni flow near the crystal and increases ½O lc again, as shown in Fig. 12. It must be noted that the
Fig. 11. Effect of gas flow rates on radial distributions of oxygen concentration at interfaces. Solid keys correspond to the results obtained without a gas guide. (a) Results for B-driven flow. (b) Results for B+Ma-driven flow.
1018 [O]lc (atoms/cm3)
Fig. 10. Radial distributions of oxygen concentration at the melt/crystal or melt/gas interfaces, mass flux and gas-phase mass-transfer coefficient with a gas guide at Q ¼ 0:5 l/min. (1) Radial distributions of [O]l (atoms/cm3) at the interfaces. (2) Radial distributions of mass flux J (mol/m2 s) of oxygen or SiO at the melt/gas interface. (3) Radial distributions of gas-phase mass-transfer coefficient, kg (m/s). Solid line: for B-driven flow. Dotted lines: for B+Ma-driven flow.
1017
10
16
0.0
B+Ma B B+Ma+Gas Guide B+Gas Guide 0.2
0.4
0.6
0.8
1.0
Q (L/min) Fig. 12. Effect of gas flow rate on the oxygen concentration at the center of the melt/crystal interface with a gas guide. Solid keys in the figure represent the results obtained without a gas guide.
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
magnitude of the Marangoni effect (i.e., temperature coefficient of surface tension, gT ) is an important parameter to determine the melt flow pattern and ultimately the oxygen concentration in this small Cz furnace.
4. Discussion The present global analyses reveal significant effects of gas-phase transport phenomena and the Marangoni effect on thermal environment, flow pattern, oxygen transport and oxygen concentration at the melt/crystal interface. However, these analyses have been conducted for a small silicon Cz furnace without rotation of the crystal and the crucible. The results may not be valid in larger Cz furnaces. The simulations will be further extended in order to include crystal and crucible rotations. Further development of the numerical code is required to enable simulations of much larger industrial Cz furnaces. In order to realize global analysis in large furnaces, we must develop a more efficient numerical code for calculating threedimensional time-dependent flow in the melt and turbulent flow in the gas phase, as well as a robust algorithm for determining time-dependent threedimensional melt/crystal interface shape.
45
flow and oxygen transport. If there were no Marangoni effect, the increase of the gas flow rate would increase the oxygen concentration in the crystal. If the Marangoni effect worked on the surface, oxygen concentration would decrease significantly by increasing the gas flow rate. These results indicate that gas flow must always be taken into account in the analysis of a Cz furnace with a heat shield. 3. Global analyses revealed that the evaporation rate of SiO is not uniform over the melt surface. In some cases, the reabsorption of SiO occurs near the crucible wall. Further simulation and development of numerical codes are needed to extend the global analyses to more realistic conditions, such as those with crucible and crystal rotations in large industrial Cz furnaces. Acknowledgements This work was supported by JSPS Research for the Future Program in the field of Atomic Scale Surface and Interface Dynamics. M.W. Li expresses his gratitude for JSPS’s support as a Postdoctoral Researcher at Kyushu University. References
5. Conclusions A global simulation code was developed to simulate the complex transport phenomena in a small silicon Cz furnace. The code could treat pseudosteady axisymmetric systems. The results revealed the effects of gas flow and the Marangoni effect on the melt flow pattern and the mass transfer characteristics of oxygen impurity. 1. The Marangoni effect (evaluated using a small temperature coefficient of surface tension, gT ¼ 7 105 N/mK) induced marked effects on the melt flow and oxygen concentration in a small Cz furnace. 2. A cylindrical heat shield, which acted as a gas guide, caused a significant decrease of heater power. It also exerted significant effect on gas
[1] N. Kobayashi, J. Crystal Growth 108 (1991) 240. [2] P. Sabhapathy, M.E. Salcudean, J. Crystal Growth 113 (1991) 164. [3] H. Hirata, K. Hoshikawa, J. Crystal Growth 125 (1992) 181. [4] T.A. Kinney, R.A. Brown, J. Crystal Growth 132 (1993) 551. [5] K. Kakimoto, K.W. Yi, M. Eguchi, J. Crystal Growth 163 (1996) 238. [6] K.W. Yi, K. Kakimoto, M. Eguchi, H. Noguchi, J. Crystal Growth 165 (1996) 358. [7] M. Watanabe, K.W. Yi, T. Hibiya, K. Kakimoto, Prog. Crystal Growth Characterization Mater. 38 (1999) 215. [8] S. Kawanishi, S. Togawa, K. Izunome, K. Terashima, S. Kimura, Jpn. J. Appl. Phys. 34 (1995) 5885. [9] I. Kanda, T. Suzuki, K. Kojima, J. Crystal Growth 166 (1996) 669. [10] N. Machida, Y. Suzuki, K. Abe, N. Ono, M. Kida, Y. Shimizu, J. Crystal Growth 186 (1998) 362. [11] D.E. Bornside, R.A. Brown, T. Fujiwara, H. Fujiwara, T. Kubo, J. Electrochem. Soc. 142 (1995) 2790.
46
M. Li et al. / Journal of Crystal Growth 234 (2002) 32–46
[12] T. Tsukada, M. Hozawa, N. Imaishi, J. Chem. Eng. Japan 27 (1994) 25. [13] H. Honda, N. Imaishi, T. Tsukada, M. Hozawa, J. Chem. Eng. Japan 25 (1992) 84. [14] A. Muhe, R. Backofen, J. Fainberg, G. Muller, E. Dornberger, E. Tomzig, W. von Ammon, J. Crystal Growth 198/199 (1999) 409.
[15] T. Carlberg, J. Electrochem. Soc. 133 (1986) 1940. [16] M. Ratto, E. Ricci, E. Arato, J. Crystal Growth 217 (2000) 233. [17] K. Nogi, private communication. [18] H. Fujii, T. Matsumoto, N. Hata, M. Kohno, K. Nogi, Metall. Mater. Trans. A 31 (2000) 1585.