Effect of Argon Flow on Oxygen and Carbon Coupled Transport in an Industrial Directional Solidification Furnace for Crystalline Silicon Ingots

Transient global simulations were carried out to investigate the effect of argon flow on oxygen and carbon coupled transport in an industrial directional solidification furnace for quasi-single crystalline silicon ingots. Global calculation of impurity transport in the argon gas and silicon melt was based on a fully coupled calculation of the thermal and flow fields. Numerical results show that the argon flow rate affects the flow intensity along the melt–gas surface, but has no significant effect on the flow patterns of silicon melt and argon gas above the melt–gas surface. It was found that the evaporation flux of SiO along the melt–gas surface decreases with the increasing argon flow rate during the solidification process. However, the net flux of oxygen atoms (SiO evaporation flux minus CO dissolution flux) away from the melt–gas surface increases with the increasing argon flow rate, leading to a decrease in the oxygen concentration in the grown ingot. The carbon concentration in the grown ingot shows an exponential decrease with the increase of the argon flow rate, owing to the fact that the dissolution flux of CO significantly decreases with the increasing argon flow rate. The numerical results agree well with the experimental measurements.


Introduction
Quasi-single crystalline silicon ingots grown by the seeded directional solidification (DS) technique are one of the main materials for solar cells, and have lower cost and a less favorable crystalline structure compared to the monocrystalline silicon. To improve silicon ingot quality, the seeded DS process still needs considerable work. The impurities in the silicon ingots significantly influence the quality of the silicon ingots and the photoelectric conversion efficiency of solar cells [1,2]. Typically, oxygen and carbon are the main impurities in the DS furnace affecting solar cell wafers. Oxygen-related defects, like thermal donors and boron-oxygen complexes, can reduce the minority carrier lifetime and cause serious light-induced degradation in solar cells [3]. Carbon-related defects, such as silicon carbide (SiC) precipitates, can lead to severe ohmic shunting, jeopardizing the electrical behavior of solar cells [4]. Therefore, it is quite necessary to control the oxygen and carbon concentrations of the DS process to produce high quality silicon ingots.
Considerable research has been carried out to control the oxygen and carbon concentrations and transport during the DS process, such as hot-zone design of furnace and optimization of operating parameters. The former method is dedicated to the configuration or the material of the local components, such as gas flow guidance device [5][6][7][8][9], graphite heaters [10,11], crucible cover [12][13][14] and heat exchanger block [15], to control the gas flow and pathway of impurity transport. The latter technique mainly focuses on modifying furnace pressure [16][17][18], gas flow rate [16,19] and growth rate [20], which is a more efficient and easier way to operate for producing high purity silicon ingots in industrial production. Liu et al. [8] pointed out that it is essential to control the argon flow in the furnace to prevent impurities from being transported to the melt-gas surface. Gao et al. [16] found that an increase in the argon flow rate can reduce both oxygen and carbon impurities in the crystals. Li et al. [19] also numerically investigated the argon flow effect on oxygen and carbon transport with a fully coupled calculation of the thermal and flow fields. However, their studies for the effects of argon flow on impurities were implemented in a small laboratory-scale DS furnace and the oxygen and carbon segregation during the DS process were neglected. Therefore, it is essential to investigate the effect of argon flow rate on the oxygen and carbon transport based on a total coupled global model to find a relative economical flow rate for industrial production.
In this study, we performed transient global simulation of coupled oxygen and carbon transport, taking into account oxygen and carbon impurities segregation during the entire solidification process. The effects of argon flow rate on the melt and gas flows, as well as the evaporation of SiO and dissolution of CO at the melt-gas surface during the growth process, were investigated. The concentrations of oxygen and carbon in the grown silicon ingot were also analyzed. The numerical results of oxygen concentration in the silicon ingot were compared with the experimental measurements acquired using Fourier transform infrared (FTIR) spectrometry. Figure 1 shows the configuration of the industrial seeded DS furnace for growing quasi-single crystalline silicon ingots. The quartz crucible has a volume of 84 × 84 × 42 cm 3 . The height of the as-grown silicon ingot is about 26.5 cm. The furnace was sealed with a water-cooled wall and inert argon gas with a pressure of 0.6 bar was used for purifying the growth environment. The furnace was equivalently simplified to be axisymmetric and subdivided into a number of block regions. Each region was discretized by structured or unstructured grids, as shown in the left part of Figure 1. For the silicon domain, a grid size of 1.5 × 1.5 mm 2 with refinement at the boundaries was adopted, which is fine enough to obtain a grid-independent solution. In-house software (CGeMoS) for modeling and simulation of crystal growth/epitaxy was used in this study. This software has been developed and used in our research [8,13,18,19] for more than 10 years. For the current grid scheme, a transient global simulation of the silicon ingot solidification process requires about 7 days on the computer with four 4.00 GHz CPUs. heaters [10,11], crucible cover [12][13][14] and heat exchanger block [15], to control the gas flow and pathway of impurity transport. The latter technique mainly focuses on modifying furnace pressure [16][17][18], gas flow rate [16,19] and growth rate [20], which is a more efficient and easier way to operate for producing high purity silicon ingots in industrial production. Liu et al. [8] pointed out that it is essential to control the argon flow in the furnace to prevent impurities from being transported to the melt-gas surface. Gao et al. [16] found that an increase in the argon flow rate can reduce both oxygen and carbon impurities in the crystals. Li et al. [19] also numerically investigated the argon flow effect on oxygen and carbon transport with a fully coupled calculation of the thermal and flow fields. However, their studies for the effects of argon flow on impurities were implemented in a small laboratory-scale DS furnace and the oxygen and carbon segregation during the DS process were neglected. Therefore, it is essential to investigate the effect of argon flow rate on the oxygen and carbon transport based on a total coupled global model to find a relative economical flow rate for industrial production.

Model Description
In this study, we performed transient global simulation of coupled oxygen and carbon transport, taking into account oxygen and carbon impurities segregation during the entire solidification process. The effects of argon flow rate on the melt and gas flows, as well as the evaporation of SiO and dissolution of CO at the melt-gas surface during the growth process, were investigated. The concentrations of oxygen and carbon in the grown silicon ingot were also analyzed. The numerical results of oxygen concentration in the silicon ingot were compared with the experimental measurements acquired using Fourier transform infrared (FTIR) spectrometry. Figure 1 shows the configuration of the industrial seeded DS furnace for growing quasi-single crystalline silicon ingots. The quartz crucible has a volume of 84 × 84 × 42 cm 3 . The height of the as-grown silicon ingot is about 26.5 cm. The furnace was sealed with a water-cooled wall and inert argon gas with a pressure of 0.6 bar was used for purifying the growth environment. The furnace was equivalently simplified to be axisymmetric and subdivided into a number of block regions. Each region was discretized by structured or unstructured grids, as shown in the left part of Figure 1. For the silicon domain, a grid size of 1.5 × 1.5 mm 2 with refinement at the boundaries was adopted, which is fine enough to obtain a grid-independent solution. In-house software (CGeMoS) for modeling and simulation of crystal growth/epitaxy was used in this study. This software has been developed and used in our research [8,13,18,19] for more than 10 years. For the current grid scheme, a transient global simulation of the silicon ingot solidification process requires about 7 days on the computer with four 4.00 GHz CPUs.  and the phase change was developed for the seeded DS furnace. The basic principle of geometry simplification from 3D to 2D axisymmetric is to keep the thermal resistance unchanged in the DS furnace. Specifically, an enthalpy formulation technique based on a fixed-grid method [21] was used to model the phase change problem of the silicon region during the DS process. After obtaining the global thermal and flow fields, coupled oxygen and carbon transport in the whole furnace based on five major chemical reactions [16,19] was calculated. In the enthalpy formulation technique, the silicon melt and crystal were treated as a unitary domain and the melt-crystal interface was not solved explicitly. Instead, a liquid fraction g l was introduced to each cell volume. Therefore, the governing equations for oxygen transport in the silicon melt and crystal are given as follows [13]:

Model Description
where ρ is the velocity, ω O is the mass fraction of oxygen atoms, u is the density, D O,m and D O,c are the molecular diffusivity for oxygen atoms in the silicon melt and crystal, respectively, and k Oseg is the segregation coefficient of oxygen in silicon. The index symbols O, m and c denote the oxygen, melt and crystal, respectively. In the current simulation, the molecular diffusivities for both oxygen and carbon in the melt and crystal are taken to be 5.0 × 10 −8 m 2 /s [22] and 5.0 × 10 −11 m 2 /s [19], respectively. Detailed governing equations and boundary conditions for the global heat transfer model have been described in our previous work [13]. The time step was set to 10 s. The governing equations for the transport of the SiO and CO species in the argon gas are as follows: where ω SiO and ω CO are the mass fractions of SiO and CO in the argon gas region, respectively, → u Ar is the argon gas flow velocity and D SiO and D CO are the molecular diffusivities of SiO and CO in the argon gas, respectively, which can be expressed as: The boundary conditions for the impurity species transport in the silicon and argon gas regions are as follows: (a) On the inside wall of the quartz crucible, the equilibrium concentrations of O atoms are expressed as: For the C atoms at the inside wall of the crucible, zero flux boundary condition is applied. (b) At the melt-gas interface, O, C, SiO and CO coexist. Equilibrium relationships for the same element between the silicon melt and argon gas, as well as the law of mass conservation, is given as follows: where c SiO and c CO are the molar concentrations of SiO and CO in the argon gas, respectively, c O and c C are the molar concentrations of O and C in the silicon melt, respectively, and R g is the universal gas constant, which is equal to 8.314 J/(mol·K).
(c) The SiO reacts with carbon at the surfaces of the hot graphite fixtures to generate CO. The reaction is assumed to be reversible and the Gibbs free energy of the reaction is expressed as follows: The equilibrium constant of the reaction is K = e −∆G/RT . The coupled boundary conditions at the graphite fixtures' surfaces can be expressed as: (d) At the argon gas inlet, the concentrations of SiO and CO are set to zero. On other solid component surfaces, zero fluxes of SiO and CO are applied. At the argon gas outlet, zero gradients of SiO and CO are used. Zero radial gradients for all the species are applied along the centerline of the fluid regions in the furnace.
The relationship between the mass fraction ω and the molar concentration c in the aforementioned equations and boundary conditions is c = ρω/M, where M is the atomic weight of the impurity species.

Effect on Heat Transfer
Argon gas is used to maintain the furnace pressure and to purify the chamber during the DS process. A proper argon gas flow rate can influence the local thermal and flow fields, which is conducive to removing oxygen and carbon impurities from the melt-gas surface for better quality ingot [23]. In this section, five cases with different flow rates (5, 10, 30, 50 and 70 L/min) at the furnace inlet were simulated to analyze their effects on the argon and melt flow during the DS process. Considering that the extreme ones, 5 and 70 L/min, are rarely used in the real industrial DS process for quasi-single crystalline silicon ingots, only the middle three cases were analyzed and presented in order to control the length of the paper. Figure 2a shows the local temperature and velocity distribution above the melt-gas surface for the argon flow rate of 30 L/min, when the solidification fraction (S.F.) of silicon melt is 50%. The injected cool argon gas flow through the centrally positioned gas inlet gives rise to a large recirculation covering the regions above the melt-gas surface and around the heaters. Figure 2b gives the velocity along the melt-gas surface for different argon flow rates. Positive values mean that the melt flows outward in the radial direction and the negative ones mean that melt flows inward. As can be found in this figure, the argon flow rate significantly influences the melt flow at the central part of the melt-gas surface. As the argon flow rate increases, argon shear stress becomes stronger and the velocity of the outward melt increases. However, the velocity shows insignificant difference at the outer part of the melt-gas surface, where the Marangoni tension dominates inward flows and the influence of argon flow is weak. The velocity along the melt-gas surface will further affect the evaporation flux of SiO and dissolution flux of CO.
L/min, are rarely used in the real industrial DS process for quasi-single crystalline silicon ingots, only the middle three cases were analyzed and presented in order to control the length of the paper. Figure 2a shows the local temperature and velocity distribution above the melt-gas surface for the argon flow rate of 30 L/min, when the solidification fraction (S.F.) of silicon melt is 50%. The injected cool argon gas flow through the centrally positioned gas inlet gives rise to a large recirculation covering the regions above the melt-gas surface and around the heaters. Figure 2b gives the velocity along the melt-gas surface for different argon flow rates. Positive values mean that the melt flows outward in the radial direction and the negative ones mean that melt flows inward. As can be found in this figure, the argon flow rate significantly influences the melt flow at the central part of the melt-gas surface. As the argon flow rate increases, argon shear stress becomes stronger and the velocity of the outward melt increases. However, the velocity shows insignificant difference at the outer part of the melt-gas surface, where the Marangoni tension dominates inward flows and the influence of argon flow is weak. The velocity along the melt-gas surface will further affect the evaporation flux of SiO and dissolution flux of CO.

Effect on Impurities Transport
Analysis of the effect of argon flow rate on heat transfer is aimed at understanding the transport mechanisms of oxygen and carbon impurities in the DS furnace. Figure 4 shows the SiO evaporation flux, CO dissolution flux and net flux of oxygen atoms along the melt-gas surface at S.F. = 50%. In Figure 4a, the result shows that the evaporation flux of SiO decreases significantly with the increasing argon flow rate. This is because the strength of the convection recirculation above the melt-gas surface becomes stronger with the increase of the argon flow rate, leading to more SiO transported back to the melt-gas surface with the recirculation which subsequently suppressed the evaporation of SiO. It is worth noting that the lowest value of the evaporation flux of SiO is located at the point of about 0.42 m along the radial position, where the outgoing and ingoing flows along the melt-gas surface meet, as shown in Figure 2b. In Figure 4b, the dissolution flux of CO also significantly decreases as the argon flow rate increases from 10 to 50 L/min, since less CO is formed and transported back to the melt-gas surface due to the decreasing of evaporation flux of SiO. To better track the transport of oxygen atoms, the net flux of oxygen atoms (SiO evaporation flux minus CO dissolution flux) along the melt surface is presented in Figure 4c. We find that the net flux of oxygen atoms along the melt-gas surface (except the peripheral region) increases with the increasing argon flow rate. The above results mean that more oxygen atoms are transported out from melt silicon and fewer carbon atoms are transported into the melt during the DS process with the increasing argon flow rate.

Effect on Impurities Transport
Analysis of the effect of argon flow rate on heat transfer is aimed at understanding the transport mechanisms of oxygen and carbon impurities in the DS furnace. Figure 4 shows the SiO evaporation flux, CO dissolution flux and net flux of oxygen atoms along the melt-gas surface at S.F. = 50%. In Figure 4a, the result shows that the evaporation flux of SiO decreases significantly with the increasing argon flow rate. This is because the strength of the convection recirculation above the melt-gas surface becomes stronger with the increase of the argon flow rate, leading to more SiO transported back to the melt-gas surface with the recirculation which subsequently suppressed the evaporation of SiO. It is worth noting that the lowest value of the evaporation flux of SiO is located at the point of about 0.42 m along the radial position, where the outgoing and ingoing flows along the melt-gas surface meet, as shown in Figure 2b. In Figure 4b, the dissolution flux of CO also significantly decreases as the argon flow rate increases from 10 to 50 L/min, since less CO is formed and transported back to the melt-gas surface due to the decreasing of evaporation flux of SiO. To better track the transport of oxygen atoms, the net flux of oxygen atoms (SiO evaporation flux minus CO dissolution flux) along the melt surface is presented in Figure 4c. We find that the net flux of oxygen atoms along the melt-gas surface (except the peripheral region) increases with the increasing argon flow rate. The above results mean that more oxygen atoms are transported out from melt silicon and fewer carbon atoms are transported into the melt during the DS process with the increasing argon flow rate.   Figure 5 shows the local distributions of SiO and CO concentrations above the melt surface (upper part) and the distributions of oxygen and carbon concentrations distributions in the silicon melt and crystal (bottom part) for different argon flow rates at S.F. = 50%. As discussed before, the vapor phase SiO reacts with the hot heater to form CO and then the reactant CO is transported back to the melt surface and dissolves into the silicon melt. Since the evaporation flux of SiO decreases with the increasing argon flow rate, both the SiO and CO concentrations decrease in the argon gas above the melt-gas surface. From the bottom part of Figure 5, we can find that the concentrations of oxygen and carbon decrease with increasing the argon flow rate. The decrease of oxygen concentration is attributed to the increasing net flux of oxygen atoms with the increasing argon flow rate, leading to more oxygen atoms taken away from the melt-gas surface, as shown in Figure  4c. The decrease of dissolution flux of CO with the increasing argon flow rate results in a decline of carbon concentration in the silicon melt and crystal. As can be seen, most of the carbon atoms were pushed into the melt due to the small segregation coefficient, which causes a relatively high carbon concentration above the melt-crystal interface and a great carbon concentration gradient along the melt-crystal interface.   Figure 5 shows the local distributions of SiO and CO concentrations above the melt surface (upper part) and the distributions of oxygen and carbon concentrations distributions in the silicon melt and crystal (bottom part) for different argon flow rates at S.F. = 50%. As discussed before, the vapor phase SiO reacts with the hot heater to form CO and then the reactant CO is transported back to the melt surface and dissolves into the silicon melt. Since the evaporation flux of SiO decreases with the increasing argon flow rate, both the SiO and CO concentrations decrease in the argon gas above the melt-gas surface. From the bottom part of Figure 5, we can find that the concentrations of oxygen and carbon decrease with increasing the argon flow rate. The decrease of oxygen concentration is attributed to the increasing net flux of oxygen atoms with the increasing argon flow rate, leading to more oxygen atoms taken away from the melt-gas surface, as shown in Figure 4c. The decrease of dissolution flux of CO with the increasing argon flow rate results in a decline of carbon concentration in the silicon melt and crystal. As can be seen, most of the carbon atoms were pushed into the melt due to the small segregation coefficient, which causes a relatively high carbon concentration above the melt-crystal interface and a great carbon concentration gradient along the melt-crystal interface. Silicon ingots were grown in the seeded DS furnace with a furnace pressure of 600 mbar and an argon flow rate of 30 L/min. The oxygen concentrations along the centerline of silicon ingot under different growth height were measured by the FTIR. Figure 6a shows the oxygen concentration along the centerline of the grown ingots. It can be found that the oxygen concentration decreases with the increasing argon flow rate due to more    1.60x10    Silicon ingots were grown in the seeded DS furnace with a furnace pressure of 600 mbar and an argon flow rate of 30 L/min. The oxygen concentrations along the centerline of silicon ingot under different growth height were measured by the FTIR. Figure 6a shows the oxygen concentration along the centerline of the grown ingots. It can be found that the oxygen concentration decreases with the increasing argon flow rate due to more oxygen atoms leaving the melt-gas surface. The calculated values of the oxygen concentration show a good agreement with the experimental measurements at the argon flow rate of 30 L/min. Figure 6b,c display the concentration and concentration differences of oxygen and carbon concentrations at Point A and Point B of the grown ingots for different argon flow rates. Point A and Point B in Figure 6 are located at 0.05 m and 0.20 m away from the ingot bottom surface along the centerline, respectively. The concentration differences of oxygen and carbon with respect to the highest argon flow rate of 70 L/min are presented. It can be found that the oxygen and carbon concentrations strongly decrease with the increasing argon flow rate when the flow rate is lower than 30 L/min, while the reduction rate becomes small as the flow rate further increases. To compromise the production cost of the argon gas and the impurity concentrations in the crystal, the argon flow rate unnecessarily increases much higher in the actual production run.  Figure 6 are located at 0.05 m and 0.20 m away from the ingot bottom surface along the centerline, respectively. The concentration differences of oxygen and carbon with respect to the highest argon flow rate of 70 L/min are presented. It can be found that the oxygen and carbon concentrations strongly decrease with the increasing argon flow rate when the flow rate is lower than 30 L/min, while the reduction rate becomes small as the flow rate further increases. To compromise the production cost of the argon gas and the impurity concentrations in the crystal, the argon flow rate unnecessarily increases much higher in the actual production run.

Conclusions
Transient global simulations of oxygen and carbon coupled transport in an industrial seeded DS furnace were performed to investigate the effects of argon flow rate on heat transfer and impurities transport. The numerical results show that the argon flow rate has significant influence on the flow intensity along the melt-gas surface. An increase in the argon flow rate can obviously reduce the evaporation of SiO and dissolution of CO along the melt-gas surface. However, the net flux of oxygen atoms away from the melt surface increases with the increasing argon flow rate. Both oxygen and carbon concentrations in

Conclusions
Transient global simulations of oxygen and carbon coupled transport in an industrial seeded DS furnace were performed to investigate the effects of argon flow rate on heat transfer and impurities transport. The numerical results show that the argon flow rate has significant influence on the flow intensity along the melt-gas surface. An increase in the argon flow rate can obviously reduce the evaporation of SiO and dissolution of CO along the melt-gas surface. However, the net flux of oxygen atoms away from the melt surface increases with the increasing argon flow rate. Both oxygen and carbon concentrations in the grown ingots exponentially decrease with the increasing argon flow rate. The numerical results of oxygen concentration along the centerline of the grown ingots agree well with the experimental measurements.

Data Availability Statement:
The data presented in this study are available within the article.