A Heat and Mass Transfer Model of a Silicon Pilot Furnace

The most common technological route for metallurgical silicon production is to feed quartz and a carbon source (e.g., coal, coke, or charcoal) into submerged-arc furnaces, which use electrodes as electrical conductors. We develop a mathematical model of a silicon furnace. A continuum approach is taken, and we derive from first principles the equations governing the time evolution of chemical concentrations, gas partial pressures, velocity, and temperature within a one-dimensional vertical section of a furnace. Numerical simulations are obtained for this model and are shown to compare favorably with experimental results obtained using silicon pilot furnaces. A rising interface is shown to exist at the base of the charge, with motion caused by the heating of the pilot furnace. We find that more reactive carbon reduces the silicon monoxide losses, while reducing the carbon content in the raw material mixture causes greater solid and liquid material to build-up in the charge region, indicative of crust formation (which can be detrimental to the silicon production process). We also comment on how the various findings could be relevant for industrial operations.


I. INTRODUCTION
SILICON is found in high quantities in the earth's crust, second only to oxygen as an element. [1] It naturally occurs in quartz and quartzite (both as silicon dioxide) and other rock, from which it can be extracted. To produce silicon, a submerged-arc furnace [1][2][3] is continually fed with carbon and quartz. [1,[4][5][6] Liquid silicon is tapped from the base of the furnace [1,4,5,7] and undergoes further processes to produce silicon with the desired purity and size distribution for the customers' needs. The hot off-gases that rise through the furnace are cooled, extracted, and filtered. Conventionally, raw materials are added in a homogenized batch at roughly 15 minutes intervals. Ideally, they should be added to balance the average rate at which material is consumed. However, in practice, a solid crust region builds up, forming a ceiling-like structure above a gas cavity. This crust is composed of a mixture of carbon, molten quartz, silicon carbide, and condensate. The presence of the crust slows the motion of solids and liquids to the lower part of the furnace, restricting their access to the heat required to facilitate the necessary reactions. The crust rises in the furnace, since material that melts and reacts away from the lower interface is replenished with additional material deposited on the top. Since the products from the lower part of the crust are tapped out of the furnace base or convected away as gas, this rising crust corresponds to a volumetric increase in the gas cavity. It is possible for the crust to be so thick and the cavity pressure so high that gas escapes rapidly out of the taphole, in a disastrous situation known as 'blowout.' [1,8] In order to counter the crust build-up, operators 'stoke' the furnace, where the crust is manually broken up by a bar from a 'stoking car' on an hourly basis. Research into silicon furnaces is thus important to improve operational efficiency and see if crust formation can be controlled.
Taking measurements in a silicon furnace is extremely difficult due to the high temperatures involved. The cost of lost production prohibits experiments on industrial furnaces. Thus, to gain insight into the furnace environment, Elkem have carried out experiments on 'Pilot Furnace Process Simulators.' These are cylinders of inner height 43 cm and inner diameter 13 cm which are filled with raw material and heated from room temperature in an induction furnace. Once the furnaces have cooled down, at the end of the experiment, epoxy is injected, so that material stays in the final state when they are cut open. Photographs of the results are shown in Figure 1. Similar experiments are described by Tangstad et al. [9] A stoichiometric model, found in References 1 and 10, captures the main reactions between furnace compounds, calculating arithmetic balances between these species to determine how much of each material is found in an inner zone at the furnace base and an outer zone at the furnace top. The model uses lumped parameters, which represent a range of physical parameters, and considers the steady state, but not how the furnace moves to this condition. The dynamic model, developed by Halvorsen et al., [4] discretizes the furnace height into several zones. Reactions take place in each zone and material and enthalpy transfers between zones at each timestep. This model is the basis for the simulator SiMod-'a Silicon Furnace Process Model.' [11] Kinetic parameters used in the dynamic model [4] and SiMod [11] are not found experimentally, but rather are chosen in order to produce sensible results that compare well with observations in industrial furnaces. Valuable background about the chemical species present in the furnace is given in Reference 12 and the relevant kinetics are discussed in Reference 13. Myrhaug et al. [6] give details about the reactions between carbon materials and silicon monoxide. Lund [14] applies existing furnace models to improve process operations. The energy balances in the silicon furnace are considered in References 15 and 16.
We develop a 1-D mathematical continuum model of a silicon pilot furnace accounting for the relevant chemistry, heat flow, and material transfer. Numerical simulations using this model are shown to agree with the results from the pilot furnace experiments, and may better inform future experimental or theoretical work involving more complicated chemical reactions or more involved geometries. Indeed, our model is modular, and certain chemical reactions can be swapped out while others added, to fit the application at hand. Furthermore, since we consider a continuum model derived from first principles, we expect that our model could be naturally extended to 2-D or full 3-D geometries.
The paper is organized as follows: In Section II, we derive the 1-D heat and mass transfer model for the silicon pilot furnace. In Section III, we perform numerical simulations in order to understand the behavior of the pilot furnace. In Section IV, we use these simulations to ascertain the influence of certain controllable parameters on the silicon production process. Concluding remarks are given in Section V.

II. MATHEMATICAL MODEL FOR A PILOT SILICON FURNACE
In this section, we develop a mathematical model which can be used to examine behavior in a silicon pilot furnace. We start by indicating which chemical species and reactions we will consider and what geometry we will use. We then formulate the framework of our model, by considering conservation of chemical species and energy, and closing the system through the constitutive relations that there are no voids, that the gas is ideal, and that the total gas pressure is constant. We then obtain the full non-dimensional model.

A. Species and Reactions
We are interested in modeling the dynamics of the following chemical species which are found in a silicon furnace: C(s)-carbon; SiC(s)-silicon carbide; SiO 2 (s)-silica (solid quartz or quartzite); SiO 2 (l)-liquid 'sticky' silica (note this can be a necessary intermediary reactant, as well as a waste slag); Si(l)-liquid silicon; CO(g)-carbon monoxide; and SiO(g)-silicon monoxide. Here, the bracketed notation (s) denotes a solid, (l) a liquid, and (g) a gas. Other compounds (such as the impurities Al, Ca, and their oxides, [1] and Ti) can be found in small quantities in industrial furnaces and may contribute to the dynamics, but are not believed to be important in the pilot furnace experiments. Fig. 1-Four pilot furnace process simulators. The furnaces were fed up to the blue horizontal line with the same 99 pct pure quartz but different reducing agents, with increasing SiO reactivity from left to right. Furnace 1 has graphitized coke which has been pre-heated to at least 2873 K (2500°C). Furnace 2 has graphitized coke and char, which has been heat treated to around 1273 K (1000°C). Furnace 3 has char only, and furnace 4 has char and charcoal, that is woodchips which have gone through a pyrolysis process at 673 K to 773 K (400°C to 500°C). The gray regions labeled droplets of silicon could also contain some slag. Some black particles in the base contain a mixture of carbon and silicon carbide, others are just carbon. Some of the most important reactions in the furnace are (see References 12 and 13) SiO 2 ðlÞ þ Si(l); ½2 We treat the chemical species as existing in one of either a joint solid and liquid bulk phase or a gas phase. We will measure chemical compositions in concentration, using the notation C X for the concentration of species X (mol/m 3 ). This is an effective concentration, since C Si measures the number of moles of Si per unit volume of the mixture of all species in the furnace.

B. Geometry
For simplicity, we model a 1-D vertical section in the furnace. Although some multi-dimensional effects are visible, such as the arched isotherms in the pilot furnace (see Figure 1) and the bridging of the crust away from the electrodes in industrial furnaces, the dominant structure is uni-dimensional. We are interested in where and how the crust occurs and this can be captured by considering height in the furnace as the spatial independent variable. By letting z denote furnace height, we take the base at z ¼ 0 and the top at z ¼ h, as illustrated in Figure 2. We allow for there to be a mixture of charge material and gas throughout the furnace, so that we do not capture the interfaces between the gas cavity and the furnace charge, or the crater and the gas cavity. Diffusion of the chemicals is neglected, meaning that transport of chemical species is assumed to be entirely due to advection of the gas phase and the solid and liquid bulk phase. It is assumed that both gases, CO and SiO, move with the same velocity v g ¼ U g k and solid and liquid components move at constant velocity v S ¼ ÀU S k, where k is a unit vector in the positive zÀdirection. Note that these are actual advection velocities, not superficial velocities, since the porosity is incorporated in the effective concentrations.

C. Conservation of Chemical Species
We define the effective mass density of species X to be p X (kg/m 3 ), which measures the mass of the species in 1 m 3 of the total mixture of all species. We define the mass density q X (kg/m 3 ) to be the mass of species X in 1 m 3 of that species alone. We assume that the mass density of solid and liquid species is constant, but the mass density of the two gases can vary. For each species, the mass density and the effective mass density can be related through a volume fraction h X (dimensionless), by p X ¼ q X h X . We note that for species X the effective mass density p X (kg/m 3 ), molar mass M X (kg/mol), and concentration C X (mol/m 3 ) are related through p X ¼ M X C X .
We consider the conservation of carbon, given by where R 1 and R 5 (both mol/m 3 s) denote the degradation of carbon in the Reactions [1] and [5], respectively. The same process can be implemented to get PDEs for the other chemical species, and we thus obtain the following system of equations:

D. No Voids Condition
For each solid and liquid species i, there is an associated volume fraction h i . There is only one volume fraction to describe the gas, h g , which will be composed of a perfect mixture of CO and SiO wherever both of these exist. If only one gas is present then the gas volume fraction is filled entirely with that gas. We make the assumption that there are no voids in the furnace, but that any arbitrary volume is filled with a combination of C, SiO 2 (s), SiO 2 (l), SiC, Si, and the gaseous CO and SiO. Thus, in terms of our volume fractions, we have that The total gas pressure in industrial furnaces has been observed to vary up to around 4 kPa [8] above atmospheric pressure (101 kPa). For simplicity, we will follow the approach used in the dynamic model of Halvorsen et al. [4] by assuming that the total pressure is constant. This constant, called P TOT , is set to 1 bar. Since a pipe allows gas to escape from the top of the pilot furnace, and the distance through which gas has to travel in the packed bed is small, this is a reasonable approximation.

F. Ideal Gas Mixture
We model the gas as being ideal, relating the gas pressure P TOT , concentration C g , volume fraction h g , gas constant R ¼ 8:314 J mol À1 K À1 , and temperature T through the equation P TOT ¼ C g RT=h g . The gas comprises a mixture of CO and SiO, so that the total gas concentration can be found as C g ¼ C CO þ C SiO . Similarly, the total gas pressure is composed solely of CO and SiO partial pressures, so that Gaskell [17] outlines Dalton's law of partial pressures, from which we can infer that We can then re-write the no voids conditions in terms of concentrations, temperature, and pressure, obtaining

G. Conservation of Energy
We now consider how the temperature is determined. In our model, we use one temperature T(z, t), which corresponds to an assumption of instantaneous heat transfer between the gas, the solid and liquid, and the graphite crucible. The temperature is influenced by heat from reactions, convection of the material in the furnace, conductive and radiative heat fluxes within the material, conduction in the graphite crucible, and external heat sources. We utilize the following equation arising from considering conservation of energy.

½17
This was derived following the approach of Reference 18 used in their model for a catalytic converter, and here, we will explain the meaning of each term in [17]. For each chemical species X , we let C kg p;X be the mass specific heat (J/kg K). We can relate this to the molar specific heat C p;X (J/mol K) through the relation C kg p;X ¼ C p;X =M X . The effective heat capacity of the solid and liquid per total unit volume is given by Although the density and specific heat of solids and liquids are assumed to be constant, r s is not constant, since the volume fractions vary. We similarly sum over the two gases to define r g :¼ h g P l q l C kg p;l , which is similarly not constant. We can write the effective heat capacities in the solid and liquid and in the gas, respectively, as The terms k f and k gr are the effective thermal conductivity of the material inside the pilot furnace and graphite crucible, respectively, and b is a radiative constant in the furnace associated with radiation between particles. Heat conducts through contact points between solid particles in the furnace and radiates across the gaseous voids between the particles. A f and A gr are the cross-sectional areas of the furnace and the graphite crucible, respectively.
The rate of heat release due to reaction j is denoted dH j , which is negative for endothermic reactions. H(z, t), the overall heat source due to the chemical reactions, is given by H :¼ P j dH j R j . The term F(z, t) denotes the external heat source going into the furnace. In the experiment, the pilot furnace starts from room temperature and is heated up to furnace temperatures by an induction heater. As can be seen in Figure 3, the bottom part of the furnace is surrounded and heated by copper coils, up to a height h T , with the top is left exposed. We take h T ¼ 0:234m from experimental data. We neglect heating into the base of the graphite crucible.
For simplicity, we assume a constant power going into the pilot furnace up to a height z ¼ h T , so that F ¼ QHðh T À zÞ, where Q is a scalar and the Heaviside function Hðh T À zÞ takes the value one for z h T and the value zero for z>h T . We let the effective power dissipated from the induction furnace into the pilot furnace be P[W]. Since Q is defined per unit volume, In the model, we have allowed convection, reaction heat release, and radiation to only act in the furnace, but all other forms of heat transfer to act in both the furnace and the graphite.

H. Reaction Kinetics
The reaction rates R j are dependent upon the local environment in the furnace, including chemical concentrations, partial pressures, and temperature. The functional dependence of the rate of each reaction was chosen after taking ideas from relevant literature [4,12,13] and the simulator SiMod. [11] Wiik [13] discusses many factors that can influence kinetic rates including the crystal structure, level of impurities, porosity, and grain size of both the quartz and the carbon used. However, we adopt the ethos of Halvorsen et al., [4] utilizing simplified kinetic rates that sufficiently capture furnace behavior. Once experimentally determined, data have been found for more sophisticated rates, these can be included without detriment to the continuum model developed in this paper.
The kinetic rates used are given in Table I. The reaction constants k j have units chosen to ensure that the overall reaction rate has units mol/m 3 s. We use concentrations to model solid and liquid reactants and a partial pressure difference from equilibrium, P SiO À P SiO;equ;j ðTÞ, for the gas reactants. Here, ðfÞ þ denotes the positive part of the function f, so that ðfÞ þ ¼ maxff; 0g. The calculations to determine the equilibrium partial pressure functions are outlined in Section II-I. Where appropriate, an Arrhenius term e ÀE j =RT is used, which is dependent on the activation energy E j , gas constant R ¼ 8:314 J/mol K, and temperature T. This is a commonly used empirical relationship to model how the rate of a chemical reaction depends upon temperature. The activation energy acts as a barrier to the reaction, the lower it is the faster the reaction will be.
In reaction R 1 , SiO gas reacts with solid carbon on the surface of particles, forming solid SiC and CO gas. Since the carbon is fed into the furnace as a raw material in the form of wood chips, coke, or otherwise, [6] some reaction models incorporate the particle size and reactivity of the carbon particles. We define this reactivity, r C , as the fraction of carbon that is available for reaction at each height and time in the furnace. One approach found in the literature to determine reactivity is to use a shrinking core model, described by References 1, 6, 19, and 20 where a spherical carbon particle reacts from the outside, with SiC product forming around a decreasing core until all the original particle has reacted. Thus, the reactivity changes with height, since carbon particles move down the furnace as they react. Andresen [12] adds further complexity by having a squared dependence on carbon reactivity, arguing that the reaction takes place in the porous volume of the carbon particles. However, for simplicity, we will treat the reactivity of carbon, r C , as a constant which takes a value between zero and one. A similar approach is taken in the stoichiometric model described in References 1 and 10 where a reactivity number, describing the fraction of carbon that will react in the upper part of the furnace, is chosen to be a constant. This approach still allows observation of the influence of more or less reactive materials on furnace dynamics.
The condensation reaction R 2 is sometimes modeled through a condensate phase Si.SiO 2 , which subsequently separates into Si(l) and SiO 2 (l), for example, in Reference 11. For simplicity, we neglect this intermediate phase and assume that SiO(g) immediately condenses into Si(l) and SiO 2 (l), which in turn can react to evaporate back into SiO(g).
Reaction R 3 represents the change of state of quartz, from a solid to a sticky liquid. We assume that solid quartz is unreactive with other substances, but that once it has melted, it may subsequently react with silicon and carbon. Alternatively, one could model SiO 2 as having solely one state. We have assumed that no SiO 2 (l) solidifies back to SiO 2 (s) and that SiO 2 (s) only melts at temperatures higher than the melting point T m ¼ 1996 K (1723 C). We model the rate of melting R 3 to be increasing with temperature and the concentration of solid quartz.

I. Equilibrium Pressure Functions
To calculate the equilibrium functions P SiO;equ;j ðTÞ , we relate the change in Gibbs free energy of a reaction, DG R [J mol À1 ], and the equilibrium constant, K equ through the reaction isotherm lnðK equ Þ ¼ À DG R RT . [17] The chemical reaction constant, K, is given as a ratio of the activity of products k þ to reactants k À , that is K ¼ k þ k À . Assuming that the activity of gases is given by their partial pressures and unitary activity of solids and liquids, we can utilize the relation P CO ¼ P TOT À P SiO to find, : ½20 The simplest means of finding DG Rj for each reaction, as outlined by Gaskell, [17] is to data fit for the constants a j , b j , and d j in DG Rj ðTÞ ¼ a j þ b j T ln T þ d j T. Using values from Reference 21, we find these constants, which are listed in Table II. We shall use these values throughout the remainder of the paper. In Figure 4, we show that [19] and [20] using data from Table II compare well with the diagram of equilibrium SiO pressures found in Reference 1.

J. Boundary and Initial Conditions
The model comprises the seven hyperbolic Eqs. [7] through [13] with reaction rates described in Section II-H, in addition to the no voids condition [16], the temperature Eq. [17], and the constant sum of partial pressures [14]. We take the constant velocity U s to be zero in our model analysis as material is not fed in from the top of the pilot furnace. We take kinetic rates k j from SiMod. [11] These were calibrated to industrial furnace data using the simplification that all activation energies E j are zero, so we follow this approach here. In addition, we have the reactivity of carbon r C , as discussed in Section II-H, as a free parameter which we can vary in our model.
To determine what boundary and initial conditions are needed for the system, we consider the characteristics of the system with direction dz dt ¼ k c in the   tÞÀplane. An explanation of how to find these is given in Reference 22. We find that there are characteristics with directions given by k c ¼ ÀU s , k c ¼ U g , and k c ¼ 1.
The direction k c ¼ ÀU s has multiplicity five, corresponding to Eqs. [7]- [11] for the evolution of the solid and liquid concentrations. Since we take U s ¼ 0, this gives k c ¼ 0 and hence z ¼ constant. Thus, we need initial conditions for each of the solid and liquid concentrations, which we take from experimental data. In the majority of the pilot furnace, there is carbon and solid quartz. However, solid quartz and silicon carbide were placed in the bottom (roughly 10 pct) of the furnace to generate the SiO furnace gas. We make the simplification that all the furnace starts with carbon and solid quartz and that there is CO and SiO initially.
Multiplying the no voids condition [16] by P TOT =RT then taking a time derivative, we can substitute in Eqs. [12] and [13] to give an equation with a spatial derivative of U g ðC CO þ C SiO Þ and a time derivative of a function of temperature and the solid and liquid concentrations. This tells us that the no voids condition [16] has characteristics with direction k c ¼ 1, corresponding to t ¼ constant. Hence, we need a boundary condition for U g ðz; tÞ. As particles touching the very bottom of the furnace remain there, we use the no-flux condition U g ð0; tÞ ¼ 0. In all practical situations, U g ðh; tÞ>0 so no boundary condition is needed at the top of the furnace.
The characteristics corresponding to Eqs. [12] and [13] are given by dz dt ¼ U g ðz; tÞ. Since we have set U g ð0; tÞ ¼ 0 then the characteristic emanating from the origin is given by z ¼ 0. Thus, only initial conditions are needed for the gas concentrations. These must satisfy the no voids condition [16], along with the initial conditions of the solid and liquid concentrations and the temperature. We are free to impose the initial partial pressure of SiO, P SiO , which we take to be 1 2 . The parabolic equation for temperature [17] requires two boundary conditions and one initial condition. Since we are heating the pilot furnace from cold, we use the initial condition Tðz; 0Þ ¼ 300 K. We assume that the furnace is suitably insulated so that there is no heat loss from conduction or radiation at the top or bottom, giving @T @z ¼ 0 at z ¼ 0; h. However, note that energy can escape from the top via convection, since gas leaves the furnace.

K. Full Non-dimensional Model
We use the non-dimensional scalings: Here, we use T bottom ¼ 2400 K (2127 C) and T top ¼ 1400 K (1127 C) as typical temperatures found in the bottom and the top of the furnace, since there is no scaling from the boundary conditions. A scale for s can be found through a balance with any of the reaction terms or the gas transport. We choose to balance with the melting of quartz reaction, so take s ¼ 1=k 3 T 0 , where we have used the shorthand T 0 :¼ T bottom À T top . This gives s ¼ 1091 seconds, which is 18 minutes 11 seconds. We take C 0 to be the initial concentration of carbon in the bulk of the pilot furnace, found to be 38, 000 mol/m 3 s. Upon non-dimensionalization, we can define the following dimensionless parameters: We drop the overbar notation and use the relations C CO þ C SiO ¼ C g , C CO ¼ C g ð1 À P SiO Þ, and C SiO ¼ C g P SiO , to write the dimensionless equations as ¼ n 1 r C C C À P SiO À P SiO;equ;1 ðTÞ Á þ n 2 À P SiO À P SiO;equ;2 ðTÞ Á þ n À2 C SiO 2 ðlÞ C Si C SiO 2 ðsÞ ðT À T Ã m Þ þ n 4 C SiC À P SiO À P SiO;equ;4 ðTÞ Á þ We have the dimensionless initial conditions (for 0 z 1) ½25 and the dimensionless boundary conditions

A. Solution Procedure
We perform numerical simulations of our model to capture the behavior in the pilot furnace. We run the simulations for 80 minutes in dimensional time, as this is the length of the pilot furnace experiments.
For our numerical simulations, we use first-order finite difference schemes, with a spatial step Dz and a timestep Dt. We first step forward in time for the solid and liquid concentrations using a Forward Euler explicit scheme. Then, we find the temperature by using a Forward Euler scheme which is implicit for diffusive derivatives but explicit for the diffusive coefficient, a þ b Ã T 3 , as well as the convective and wall source terms. The reaction heat release is implicit in solid and liquid concentrations but explicit in temperature and pressure. We use a Godunov-type method with upwinding for the convective fluxes. In the time partial derivative, we neglect the heat capacity in the gas r g , which is 10 À4 times the size of the other heat capacities r s and r Ã gr . Using the temperature, we then find the concentration of gas using the no voids condition [21]. We next update for the gas velocity and partial pressure of SiO using a Godunov-type upwinding method. In order for stability, we choose Dt based on Dz to obey the Courant-Friedrichs-Lewy (CFL) condition. The scheme implemented is grid scalable and tests show a slightly better than first-order rate of convergence.

B. Model Behavior
Typical dimensionless parameters are given in Table III. These were calculated from dimensional parameters found in References 1, 11, 21, 23 through 29; for brevity, we do not list all of these values here as they are never used individually. Unless stated otherwise, we use the value r C ¼ 0:2 for the reactivity of carbon in simulations.
In the first 40 minutes, the furnace heats up and no substantial reactions occur. Then, in the second 40 minutes, the pilot furnace is hot enough for typical industrial furnace behavior to be observed. A separation emerges between a lower region, where SiO 2 (s) melts and carbon is degraded, and an upper region, where the SiO gas flowing up from the bottom condenses. The interface between these regions moves upwards as the temperature rises throughout the furnace and more of the quartz becomes hot enough to melt. Plots of temperature and the total solid and liquid volume fraction h S are given in Figure 5.
In Figure 6, we compare our numerical simulation with one of the pilot furnace experiments, displaying results in terms of volume fractions. In the top region, z ! h Ã T % 0:54, we have unmelted SiO 2 (s). Carbon has reacted away in part and there is deposition of SiC, Si, and SiO 2 (l). The SiO 2 (l) comes from condensation of SiO, rather than melting of SiO 2 (s).
In the bottom of the furnace, z h Ã T % 0:54, SiO 2 (s) has melted and there has been degradation of carbon, through reaction R 5 with the SiO 2 (l). There is some build-up of SiO 2 (l) and SiC, as well as production of Si. Around the highest position of heating z ¼ h Ã T % 0:54, in approximately the region 0:45 z 0:6, the volume fractions of Si and SiO 2 (l) seem to interchange between their behavior in the upper furnace and lower furnace. This crossover region corresponds to the temperature around 2000 K (1727 C), where both quartz first starts to melt and the SiO gas begins to condensate.
Our 1-D model captures a distinct upper charge region, and shows the interface of this region moving upwards as temperature increases. To replicate the bridging of the crust, a 2-D model would be necessary. In the experiment, there is a separation between the gas cavity and the crater comprising of solid and liquid. This is not replicated in our model, as we treat solid and liquid as being stationary. In the future, we could model the dripping of material from the base of the furnace charge into the crater, as well as the flowing of the highly inviscid Si. In Figure 7, we plot the solid and liquid concentrations and the reaction rates showing the dynamics in the pilot furnace at 40, 60 minutess, and the end time of 80 minutes. The temperature increases with time throughout the furnace in a monotone profile which is hotter at the base and cooler at the top. There is a separation of regions in the furnace with the amount of solid and liquid concentrations. SiO 2 (s) melts in the bottom and forms SiO 2 (l), which reacts with carbon to produce CO and SiO gas. This can be seen through R 3 , the melting of quartz, and R 5 , the carbon and SiO 2 (l) reaction, which dominate in the base. The height of this base region increases with the temperature. In the top of the furnace, SiO gas condenses to form Si and SiO 2 (l). The SiO 2 (l) then reacts with carbon to produce gas, as occurs in the furnace base. Some of the SiO gas reacts with carbon to produce CO gas and SiC.

IV. INFLUENCE OF MODEL PARAMETERS ON SILICON PRODUCTION
In the silicon process, and the pilot furnace experiments, various inputs can be changed, including the reactivity of the carbon particles used, the heating applied to the furnace, and the ratio of carbon to quartz used as raw materials. These are varied in our model to see the influence on furnace dynamics. Similarly, although the kinetic rates cannot be changed in reality, since they are not known accurately they can also be adjusted in our simulations. Of particular interest is the condensation reaction R 2 , so we vary the kinetic rate n 2 .
We want systematic ways of measuring how furnace behavior changes with varying inputs, and so use the metrics of silicon yield, SiO losses, and solid build-up, described below.
We define the silicon yield (percent) as follows: Silicon Yield (percentÞ ¼ Total Final Silicon Total Input SiO 2 Â 100; ½27 where we measure in terms of concentrations. This calculates the amount of silicon that is created from the raw material quartz. The more silicon that is created, the better. Due to the uncertainty of kinetic parameters and purity of quartz used, these calculated values are not intended to be used to a high level of accuracy, but are useful for comparison between different model inputs. The values will be lower in the pilot furnace experiment than in industrial furnaces, because the experiment only runs for 80 minutes and quartz in the upper furnace remains unmelted and unreactive, so results should not be directly transferred to the operational setting. We define the SiO losses as SiO Losses (percent) ¼ Total SiO convected away Total Input SiO 2 Â 100;

½28
where again we measure in terms of concentrations. This determines how much of the silicon put into the furnace (as SiO 2 (s)) is lost from the process through SiO gas escaping out of the furnace top, but does not include losses in the tapping process or during post treatment. It is an important measure, capturing the efficiency of production, since the majority of silicon that is not lost from industrial furnaces through SiO convection will eventually be tapped as liquid silicon. We want to have a measure for crust formation and so, as a proxy, consider the build-up of solids and liquids. We define the solid build-up as where h S is the total solid and liquid volume fraction and we are thus finding the biggest increase in h S from the initial condition. The higher this value is, the more clogging up of solids and liquid we see in the upper part of the furnace. We run simulations with varying inputs and calculate the output metrics described. In Figure 8, we show plots for varying the reactivity of carbon and applying three different levels of effective heating.
For the higher heating powers 15 and 20 kW, there is a trade-off between silicon yield and SiO losses. Higher heating brings higher yield, but also higher losses, as more SiO escapes from the system. It is difficult to apply these results to industrial furnaces, where there is feeding of raw material and the heat source is from electrodes, not an induction heater. However, if the power input is too high, it is possible that much of the silicon will be lost from the process as gas. For all levels of heating in the pilot furnace, SiO losses decrease as the carbon reactivity is increased, since reaction R 1 between SiO and carbon is quicker, so less SiO escapes the furnace, and more SiC is produced. For a heating of 10 and 15 kW, the solid build-up is approximately constant for varying reactivity of carbon. Higher reactivity causes greater SiC production, which takes up more space than the carbon it replaces (the volume taken by two moles of carbon is 2 Â M C =q C ¼ 1:06 Â 10 À5 m 3 , but the volume taken by one mole of SiC is M SiC =q SiC ¼ 1:25 Â 10 À5 m 3 ). However, an increase in R 1 produces more CO(g), lowering the partial pressure P SiO . This lower P SiO decreases the rate of R 2 , reducing the deposition of SiO 2 (l) and Si, and balancing the build-up of SiC.
In Figure 9, we show plots for varying the mixture between carbon and solid quartz used at the start and applying four different values for r C . The mixture ratio is the initial value of C C =C SiO 2 ðsÞ . Note that when we increase the carbon, we also decrease the quartz, in order to ensure that the initial values of h S are the same. Thus, the silicon yield and SiO losses are measured with respect to the initial amount of quartz, which is different in each simulation. We see that more reactive carbon leads to a lower silicon yield and lower SiO losses, since more SiC is produced through the reaction R 1 and less SiO can escape from the top of the furnace. However, silicon yield increases with reactivity for effective heating of 15 and 20 kW (see Figure 8). Also, in an industrial furnace, most of the silicon atoms remaining in the process will eventually be tapped as liquid silicon, so SiO losses could be a more useful metric than silicon yield.
We see that as the concentration of carbon in the initial mixture increases (and correspondingly the amount of solid quartz decreases), there is a higher silicon yield and lower level of SiO losses, due to a faster rate of R 1 . As with Figure 8, this faster R 1 decreases the rate R 2 , due to lower P SiO , causing less SiO 2 (l) and Si to be deposited. However, in this case, this effect dominates over the build-up of SiC, and hence the solid build-up decreases with more carbon in the initial mixture.
In Figure 10, we see the effect of changing the strength of the condensation reaction R 2 in our simulations, by altering the kinetic rate parameter n 2 . We find that the silicon yield increases, SiO losses decrease, and the solid build-up increases with a higher value of n 2 . This is expected, since faster R 2 means a greater rate of condensation of SiO gas to SiO 2 (l) and Si.

V. CONCLUSIONS
We have developed a mathematical model for the chemistry, heat flow, and material transfer in a silicon furnace, and we have applied this model to better understand silicon pilot furnace experiments. For higher process yield and energy efficiency to be possible in industrial silicon production, it is important to understand crust formation within the furnace. Numerical simulations replicated the reacting away of the lower part of the furnace, which corresponds to the creation of a gas cavity and crater region. Comparisons with experiments show a good fit for the position of the interface between the charge and the lower regions, without the need to fit any model parameters artificially. The location of this interface is largely dictated by temperature, with our results suggesting that increasing temperature results in an increase in the size of the gas cavity and a rising of the crust within the furnace. Key to modeling the furnace dynamics is understanding the gaseous SiO molecules convecting up the furnace. Due to the partial pressures and temperatures found in the upper region of the furnace, these molecules can either condense to form Si and SiO 2 (l), react with carbon to form SiC and CO, or convect out of the furnace (resulting in silicon loss). It is primarily the balance of these three processes which determines the productivity of the silicon furnace. In furnace operation, the SiO losses and crust formation can be reduced simultaneously. Increasing the mixture ratio and the reactivity of the carbon particles decreases the crust formation, seen through the proxy of solid build-up and the amount of SiO convecting out of the furnace. However, using highly reactive carbon particles is costly, so the financial benefits of reduced SiO losses should be considered when determining what constitutes optimal production.
Our simulations suggest that material builds up in the upper region of the furnace due to the production of SiC and the condensation of SiO gas. Having more carbon in the initial mixture plays a larger role in the reduction of the solid build-up than does using more reactive carbon particles. Yet, using too much carbon and too little quartz may reduce the overall amount of silicon produced, even if the silicon yield and SiO losses are favorable, as these are measured relative to the amount of initial quartz. On the other hand, using more reactive carbon particles decreases both silicon yield and SiO losses for the effective heating used in the experiment, as more SiC is produced. Since SiO losses is a more helpful metric for industrial operation, the finding that these diminish as reactivity of carbon increases confirms the intuition presented in Reference 1.