Numerical modelling of surface waves generated by low frequency electromagnetic field for silicon refinement process

One of the most perspective methods to produce SoG-Si is refinement via metallurgical route. The most critical part of this route is refinement from boron and phosphorus, therefore, approach under development will address this problem. An approach of creating surface waves on silicon melt’s surface is proposed in order to enlarge its area and accelerate removal of boron via chemical reactions and evaporation of phosphorus. A two dimensional numerical model is created which include coupling of electromagnetic and fluid dynamic simulations with free surface dynamics. First results show behaviour similar to experimental results from literature.


Introduction
The exhaustion of fossil fuel resources, global warming, and increased energy consumption all make finding alternative energy sources important. In terms of environmental impact, solar energy is one of the most perspective and attractive sources. However, development of solar power is limited by the solar power costs in comparison to other power sources. As shown in [1], by 2020 the energy cost for solar power is still expected to be higher than for hydro, onshore wind, geothermal, nuclear and even biomass energy sources. Reduction of solar cell production costs is necessary for this technology to develop, and one solution is to reduce the cost of raw materials like polysilicon.
In 2013, production of polysilicon constituted approx. 23% of the total solar cell costs [2]. This number can be reduced by lowering energy consumption during solar grade silicon (SoG-Si) production. One of the most efficient ways to produce SoG-Si is the metallurgical route of refining. It consumes about 30 kWh/kgSi (kilowatt-hours per kilogram of silicon) of energy. In comparison, chemical route such as the Siemens process, which is the most widely used refining process nowadays, consumes more than 70 kWh/kgSi [3]. In 2013, more than 200 000 tons of polysilicon was produced worldwide, which could result in more than an 8 TWh reduction if all of the silicon was produced via the metallurgical route rather than with the Siemens process.
Currently, the polysilicon market is recovering from crisis and new technologies that ensure reduction of production costs can give a significant advantage in a competitive SoG-Si market. For this purpose sophisticated simulation tools are required, due to complexity of underlying phenomena, and multiphysical nature of this process, which includes electromagnetics, thermal physics and fluid dynamics. Free surface phenomena is of particular interest, because of its complexity and huge computational requirements due to necessity of coupling between different simulation tools or packages. There are several papers on free surface dynamics in electromagnetically driven flows published in recent years. In [4] it is shown that volume of fluid (VOF) method gives good agreement with experimental results. Drawback of use of VOF for problems for very dynamic surfaces is large computational time. On the other hand, when surface deflection is not large, VOF allow calculation of free surface dynamics in reasonable time [5]. In all mentioned works, it is common assumption that EM frequency is too large to influence surface dynamics at its frequency. However, in problem under consideration in this work, it is not the case. In this work low frequency EM field is used and time dependence of EM variables has to be taken into account. The main task of this work is creation of numerical tools for simulations of silicon refinement from boron and phosphorus with use of advanced electromagnetic impact.

Potential for EM technologies in silicon refinement
In many investigations phosphorus removal is found to be the most challenging task, since the purity threshold (< 1 ppm) was not reached [6]. Although, oxidation and slag refinement is simple process for removing phosphorus from silicon, it is insufficient to reach low purity threshold. Mostly used approach for phosphorus removal is evacuation, it is dependent on free surface area and bulk melt stirring. Refinement enhancement can be achieved by superposing DC and AC magnetic fields [7]. Furthermore, this impact creates also mixing effect in the bulk melt, which is also important for boron removal. In comparison, classical induction stirring can enhance only bulk melt mixing, leaving surface area size unchanged.
Boron removal by oxidizing bubbles was already shown in [8] where significant splashing of silicon was observed, as a result blowing of oxidizing gas on free surface of the silicon was preferred instead. Also, no attempt was made to adapt electromagnetic stirring for bubble flow control, which was shown to be possible in [9]. Furthermore, high velocity shear in melt can lead to gas bubble breakup thus enlarging gas-liquid interface area and increasing boron removal rate.
Other possibility for boron removal could be blowing of oxidizing gas on free surface, which is already well known [10], and with additional creation of surface waves by means of electromagnetic impact thus creating capillary waves and increasing surface area at least twice. Previously mentioned electromagnetic impact technology (superposition of DC and AC fields) can be used for boron evaporation, because it enhances surface area and stirring.
Highest boron refining rates can be achieved with plasma refining, but it has high energy demand. In reactive gas refining, concentration of boron changes according to first order rate law (1) Here CB -boron concentration with indices after comma indicating time, A is the area of interfacial surface between oxidizing gas and silicon melt, V is the volume of silicon melt, t is time, kB is boron mass transfer coefficient at the surface, measured in μm/s. kB is determined experimentally and is dependent on oxidizing gas supply rate, type of oxidizing gas and melt temperature. Creation of surface waves would increase A/V ratio, thus speeding up boron removal. Table 1 summarizes kB from literature with different oxidizing gases, gas flow rates and temperatures.
Data from Table 1 shows that the most important factor for boron removal is oxidizing gas and its flow rate.
According to [13], the same equation (1) is applicable for phosphorus removal, literature data is also shown in Table 1. For phosphorus, temperature is a more important factor. Another important factor is pressure. According to [14], the phosphorus mass transfer coefficient on the surface of the silicon melt strongly depends on pressure. Above 1 Pa pressure, the transfer rate of phosphorus declines very rapidly. On the other hand, phosphorus mass transfer coefficient reaches asymptotical value below 0.1 Pa and therefore further decrease in pressure would not be beneficial.
For removal of both boron and phosphorus removal, surface waves can be used to enlarge surface area and increase refinement rate. However, it is not clear how large the surface enlargement will be and 3

1234567890''""
VIII International Scientific Colloquium on Modelling for Materials Processing IOP Publishing IOP Conf. Series: Materials Science and Engineering 355 (2018) 012020 doi:10.1088/1757-899X/355/1/012020 how it will influence boron diffusion in gas boundary layer. For this purpose, the first step is creation of a numerical model which can simulate melt's surface waves induced by electromagnetic impact.

Simulation model
A 2D axisymmetric problem of mercury with free surface is considered. The mercury tank is surrounded by inductor, which creates alternating magnetic field, see Figure 1 (a) for illustration. A similar device and experiment was described in [15]. This mercury experiment is only used for model validation purposes.
In For numerical simulations we use two open-source software packages -the finite element multiphysics simulation code Elmer and computational fluid dynamics code OpenFOAM. Efficient parallel coupling between both codes is done using EOF-Library [16].
Free surface is computed using Volume of Fluid (VOF) method, which is implemented in interFoam solver in OpenFOAM. Euler scheme is used for time discretization, linear schemes for spatial discretization of all other terms in momentum equations. Equations for turbulent properties k and ε are discretized using upwind scheme.
Computational mesh for OpenFOAM has 51k cells, mesh for Elmer has 46k elements. We used very fine mesh with cell size 0.085x0.35mm in the region where surface can be present during simulation. Computation scheme is shown in Figure 1 (b). OpenFOAM package solves fluid dynamics with free surface using volume of fluid method, then computes electrical conductivity of gas-and-melt mixture and sends it to Elmer package. Elmer solves the complex problem for the time-harmonic electromagnetic (EM) field and sends back to OpenFOAM the complex components of electrical current density and magnetic field density.
Every time the EM field is updated, the amplitudes are computed (2) and the corresponding phases for current and magnetic field densities, where i represents x, y and z components (radial axial and azimuthal respectively), and atan2 is the arctangent function with two arguments that gives the angle in radians between the x-axis and the specified point. Then, at every OpenFOAM timestep the Lorentz force is updated where ω is angular frequency and Fcorr is a force correction due to interaction of velocity and magnetic field: An exponential multiplier in equation (4) provides a smooth increase of Lorentz force during the first few periods and reduces unwanted splashes at the beginning of simulation.

Results
Turbulent flow (Re ≈2.5•10 4 ) inside mercury vessel has two toroidal vortices similar to a classical induction crucible. Steel ring near the bottom of the crucible (see Figure 1 (a)) significantly changes the size ratio of two vortices -the upper vortex becomes larger and more dominant due to field screening near bottom of crucible. Figure 2 (a) shows instantaneous velocity field in the melt. The velocity vectors in this figure that have also vertical component only indicate that melt surface is moving and there is not net melt mass moving through the interface.
Slow change of the EM force in the melt causes creation of surface waves near the outer wall, and subsequent propagation of waves radially inward. Waves have very small amplitude (below 0.5 mm) and they are hardly distinguishable in Figure 2(a) Frequency analysis of local melt height in three points (x = 0 on axis, x = 35 mm, and x = 70 mm near wall) is shown in Figure 2 (b). There is one very strong peak at 28 Hz that is connected with the frequency of the EM force, which is double the EM field frequency. Analysis of surface data shows that there are 10 full waves on the surface, which gives wavelength λ ≈ 7 mm, same result as in [15]. The time dependence for the surface shape of melt during one 14 Hz period can be seen in Figure 3 (a). It was observed that the surface has regular structure with waves traveling radially inward. These data also allow estimation of wave phase speed by means autocorrelation of point data. Obtained result is c ≈ 19.5 cm/s, and it corresponds well to λ/T=19.6 cm/s.
Another comparison was done with experimental velocity measurement data from [15]. Figure 3 (b) shows radial velocity component plot on vertical line r = 49 mm. There is a good agreement between experimental and simulation data with ring. Initially simulations were done without steel ring, which was not mentioned in [15], but was found in other work of the same author [17]. This result shows how strongly steel ring changes the flow structure. There is uncertainty about wave height in experimental data [15] and therefore validation of model requires further step. Furthermore, obtained results show good agreement, but obtained wave heights are too small for practical applications of silicon refinement. For this purpose, the magnetic field on crucible axis was increased from 0.108 T to 0.257 T.
For this case, wave structure differ significantly ( Figure 4). Wave amplitude is much larger (up to 2 mm), but there is no regular structure of waves moving toward centre. Waves here can be separated in three regions.
Most outer region, between radius values 45 and 70 mm, is characterized by waves traveling radially inward. However, they differ from small wave case, because here one wave is created for every two force pulses. In small wave case one waves was created each force pulse.
In middle region (15-45 mm) waves are travelling radially toward the centre initially and when waves reach the centre, waves travelling away from centre appear. After two seconds, first standing waves start to appear when oppositely travelling waves interfere.
In innermost region spatial waves are almost not distinguishable, but whole region oscillates in time.
It is visible that for high magnetic field case, wave height strongly decrease in radial directionthere is approx. 2 mm amplitude near outer wall and about 5 times smaller amplitude near the centre of crucible.
This finding with higher amplitude waves correspond to experimental findings in [15] where at high magnetic field, structure of waves becomes non-axisymmetric. For this reason, extension of this work with 3D simulations is planned.

Conclusions and outlook
Developed numerical model showed that waves are formed on melt's surface with low frequency EM field excitation. Results are in good agreement with literature data. However, these waves are too small to be used in practice.
It was shown that waves with larger amplitude differ significantly in structure and therefore require more thorough study and experimental validation. There are very few references on numerical modelling of such problems and therefore these promising results are important step towards better understanding and capturing of EM induced waves on liquid metal surface.