(1) Institute for Electrothermal Processes, University of Hannover
Wilhelm-Busch-Str. 4, D-30167 Hannover, Germany

(2) Laboratory for Mathematical Modelling of Environmental and Technological Processes, University of Latvia
Zellu str. 8, LV-1002 Riga, Latvia


Numerical modeling of the concentration and temperature distribution in axial symmetrical systems with several recirculated flow eddies, which is based on various 2D stationary k- e models and commercial codes, e.g. ANSYS and FLUENT, leads to results, which are significantly different from experimental data. Therefore additional user-defined subroutines were included in the commercial program code to improve the turbulent heat and mass transfer in the zone between the recirculated flow eddies. In addition transient 3D calculations were performed to investigate scientifically the flow dynamics.

1. Introduction

In many electrometallurgical applications, for example, induction crucible furnaces (ICF), an intensive turbulent melt flow up to 1 m/s with several typical eddies of the main flow is formed under the influence of the electromagnetic forces. Measurements of the turbulent flow velocities in the experimental ICF show large intensity of the velocity pulsations especially in the region between the macroscopic toroidal flow eddies [1]. These oscillations are oriented mainly perpendicular to the time-averaged flow and therefore the heat and mass exchange between these two regions is caused by this pulsation. Experimental results show, that there are intensive low-frequency oscillations with a time period from 5 up to 20 sec, depending on the crucible geometry and the time-averaged velocity of the recirculated flow, in the zone of the eddies interaction [2]. Traditional numerical calculations of the flow and the heat and mass transfer processes, based on commercial codes using various stationary 2D k- e models show that these models do not take into account the low-frequency oscillations of the flow [3,4]. Therefore the calculated temperature and concentration distribution is essentially different from experimental results, although the predicted averaged flow and typical values of velocities are in a good agreement with the measurements. Taking into account the technological importance of the numerical modeling of these widespread industrial processes, two ways for solving this problem were realized:

2. Analysis of experimental data

The experimental investigations of the velocity and temperature distribution were carried out using a crucible inductor furnace (CIF), which consists of a double-walled cylindrical vessel, where an ICF is flanged below. Wood’s metal, which has a melting point of 72 °C was used as a model melt. The local flow velocities within the melt were measured using a well-proved measurement system based on a magnetic sensor [5]. Distances between the local measurement points were 0.3 – 1.5 cm in radial and 3 cm in axial direction (4x4 cm in upper part). Figure 1 shows the characteristic time-averaged velocity distribution with two toroidal eddies of the recirculated flow formed under the influence of electromagnetic forces when operating frequency is about 500 Hz. Maximum velocity is approximately 45 cm/s, therefore the flow is highly turbulent (Re »105 ). In the upper part of the installation (furnace vessel) the electromagnetic field intensity is nearly zero and the flow is mainly driven by buoyancy forces in that region, resulting in small velocity values (under 5cm/s). Comparing the calculated and measured velocities it must be considered that maximal values are located very close to the crucible wall, where they can’t be precisely measured with the sensor. That’s why, it is reasonable to take characteristic values on the symmetry axis.

Induced heat sources are distributed almost equally in the lower and upper eddies. The averaged temperature difference measured between the lower and upper flow eddy is only about 5 K, and this shows, that the heat exchange is well developed in the interaction region.

The melt flow in induction furnaces was many times simulated numerically with the help of the k- e and similar two-parameter turbulence models, which are widely used for various engineering applications. These calculations usually produce good results for the averaged velocities, but the calculated temperature distribution very often is different from the experimental results [3,4,8]. For example, we present our result for 2D axial symmetric system, where flow calculations were based on RNG modification of the previously mentioned k- e turbulence model, which is implemented in commercial CFD software FLUENT. In Figure_2 is shown, that calculated and measured velocity distributions are in a good agreement. However the temperature fields have principal differences, and the main distinguishing feature is very large temperature gradient between the lower and upper toroidal flow eddy when calculating numerically ( Fig. 5a ). Within the k- e turbulence model turbulent viscosity is calculated in the following manner:


where k is the turbulent kinetic energy, e is the dissipation rate of k and C m is the model constant. In Figure 3 is shown, that the turbulent viscosity is one order of magnitude higher than the molecular viscosity in the region of eddies interaction. In addition, the axial velocity is nearly zero in the same region. Therefore, both, the diffusive and convective transfer mechanisms do not provide sufficient heat exchange between the lower and upper toroidal flow eddies. As result, the simulation predicted larger temperature difference than was measured in the experiment ( Fig. 5a ).

The k- e model calculates the turbulent thermal conductivity with following expression [7]:

where Prt is the turbulent Prandtl number and cp is the specific heat. The calculations show that maximum of turbulent viscosity is located in the centers of the eddies ( Fig. 3 ), therefore the maximal thermal conductivity is located there too. Even if we reduce the turbulent Prandtl number, for example Prt=0.1 instead of standard value 0.85, it is not possible to improve the turbulent heat transfer between eddies. That’s why, k- e turbulence model can’t correctly describe transfer processes in this case.

Trying to understand the reason why numerical and experimental results are so different, transient measurement data for the velocity components were analyzed. These velocities were measured with a scanning rate of 20 s -1 during 50 seconds. It was found out, that the flow intensively oscillates and these oscillations have low-frequency components. Figure 4 shows the local axial velocity component’s dependence on time in a point, which is located between the eddies near the crucible wall. Such velocity’s behavior is characteristic for all other points in this region. The amplitude of the oscillations is about 30 cm/s and this value is in the same order of magnitude as the maximum time-averaged velocity, so the influence of this phenomenon on the heat and mass exchange processes has to be taken into account. Auto-correlation analysis shows that the long-time period in this case is about 3.5 seconds. We can imagine that the heat transfer is realized by the melt moving between the lower and upper region. The melt from the lower region flows into the upper region to be cooled down there, then the melt moves back to be heated up again. The steady-state simulations, of course, can not take into account such transfer process because of its long-time character. The amount of heat, which can be transferred by the described mechanism, can be estimated by the melt movement through the virtual border between eddies in a unit of time. In this way the heat flux F can be described by the moving mass of the fluid, the average temperature difference DT between the lower and upper region and the specific heat cp:


The mass flow can be calculated considering the intensity of the axial velocity oscillations, which changes in radial direction, where near the symmetry axis the oscillations are significant smaller. The result of this estimation is a total heat flux of about 9.5 kW between the lower and upper flow eddy. This estimation is in good agreement with the numerically calculated and experimentally measured induced power in the lower eddy of about 10 kW. Taken into account, that heat losses through the crucible wall are insignificant the described transfer mechanism plays the main role in heat exchange processes between the averaged flow regions.

An effective thermal conductivity can be derived from the calculated heat flux:


The effective thermal conductivity can be described by the sum of the turbulent conductivity l t, the molecular conductivity land the conductivity l lf provided by low-frequency oscillations:


It is noticeable, that l lf is one or two orders of magnitude higher than the calculated with k- e model turbulent and molecular conductivity in the examined region.

3. Numerical modeling with additional thermal conductivity .

The development of an extension for the standard k- e turbulence model, which would take into account the macroscopic heat transfer mechanism in stationary calculations of the fluid flow can be realized in several ways. One of the approaches is to take into account the effect of low-frequency oscillations in the turbulent kinetic energy generation term in the k- e model [4], This modification has significant influence on the turbulent kinetic energy distribution. But these long-period oscillations are not actually turbulent, and adding their kinetic energy to the generation term results in increasing of the turbulent viscosity. It results in changes in the averaged flow and enlargement of the zone between the flow eddies. This doesn’t lead to the improved heat transfer.

Another way is to increase the effective thermal conductivity accordingly to the values achieved from the experimental data. This was realized by a user-defined subroutine implemented in the commercial program FLUENT, which dynamically, during each iteration cycle, recalculates thermal conductivity as described below.

Considering the large amplitudes of the velocity oscillations, the distribution of the effective thermal conductivity can be assumed almost homogeneous in the region of eddies’ interaction when calculating stationary. The calculated turbulent conductivity in the centers of eddies is taken as the base values. To describe the distribution of this flow parameter in the intensively oscillating flow, the linear interpolation of these values in the axial direction between eddies’ centers is used. Also, the turbulent Prandtl number was changed to appreciate the influence of the convective compared with the diffusive heat transfer. A new value of 0.16 instead of 0.85 in the standard model was obtained for this parameter from numerical studies.

This extension was tested for different sets of experimental data (various currents in the inductor and frequencies) and the calculated temperature interval is in good agreement with the experiment ( Fig. 5b ). Also, this method can be used to calculate the mass transfer processes in induction furnaces, for example, in situations when alloying compounds are added on the surface of the melt.

4. Numerical modeling of the eddy dynamics.

The numerical investigation of the velocity oscillations requires 3D transient analysis, because of three-dimensional character of this phenomenon. The crucible inductor furnace (CIF) is not suitable for studies of such kind of modeling, because in the furnace vessel the flow is mainly driven by buoyancy forces and this situation can significantly increase the convergence time because calculations of flows with natural convection usually converge very slow. In comparison with the CIF the standard ICF has a simple geometry suitable for first studies of 3D numerical modeling. Experiments show that this crucible has practically homogeneous distribution of the temperature, so the thermal convection can be neglected in order to save computational resources. A large amount of experimental data is available for this furnace [5]. The measurements show, that in the melt of the ICF also exist low-frequency velocity oscillations with long time periods between 9 and 12 sec depending on inductor current I as shown in the Figures 7 and 8 . This dependence is similar with the expression for period of inertial waves in rotating fluid when the wave front is moving in the rotation axis direction [6]:

,     (4)

where W is an angular velocity. This angular velocity in particular point can be calculated as

,     (5)

where vch is the characteristic flow velocity and rchis the distance from the eddy centers. Considering that vch ~ I, we can conclude T ~ 1/ I ( Fig. 7 ). It must be mentioned that auto-correlation analysis of the oscillations in other regions of the crucible shows different periods due to the influence of the crucible wall and the free surface of the melt.

The fact, that such phenomenon can be predicted using basic transient three-dimensional hydrodynamic equations, gave the premises, that similar oscillations could be obtained using standard commercial 3D CFD software, e.g. FLUENT. This program performs hydrodynamic calculations basing on the common fluid flow equations, which are solved in an iterative manner with the help of a numerical algorithm. It must be considered that results depend on numerical realization in a certain extent.

Calculations were based on the standard k- e turbulence model. Mesh contained about 230,000 tetragonal elements inside the crucible. The simulation converged to the almost stationary state with symmetric recirculated flow eddies. Considering, that the time averaging of the flow, during iterations at each time step, is able to damp flow perturbations, the time-step interval was reduced from 0.05 sec to 0.02 sec, but the flow remained stationary without oscillations.

In the real crucible long period oscillations can arise due to the flow instabilities, which are caused by asymmetrical inductor and other perturbations. Considering this, there were attempts to swing the flow by raising the relaxation factors used by the solver, but often this led to the divergent solution. One of the reasons which may damp the oscillations, is the presence of relative high values of the turbulent viscosity (up to 1000 times higher than molecular viscosity) in the centers of the eddies. By the way, such distribution has no physical explanation.

Another possible approach to transient numerical simulation is LES (Large Eddy Simulation) turbulence model, which is considered to be some kind of compromise between the k- e model with relative low mesh quality requirements and DNS (Direct numerical simulation) method based on non-averaged Navier-Stokes equations. Mesh with more than 3.5 millions elements was generated to try out this scheme. Transient calculations with time a step of one second produced quasi-stationary state of the averaged flow ( Fig. 9 ), which is in good agreement with the experimental data. But using LES the distribution of eddy viscosity essentially differs from results received by the use of k- e model ( Fig. 10 ). Large eddies are resolved directly, therefore subgrid turbulent viscosity characterize dissipation of eddies, which are smaller than grid element size. As result, maximal values are one order of magnitude smaller than in case of k- e model, but now these maximums are located precisely in the zone of macroscopic toroidal eddies interaction as shown in the experiment. Continuing these calculations with small time step of 10 ms, flow oscillations were obtained, which have similar parameters with derived from experimental data. Transient three-dimensional LES numerical results will be presented in our future publications.

5. Conclusions

The calculations, which were performed using the extended stationary 2D k- e turbulence model taken into account the additional heat transfer mechanism by the long period oscillations, lead to results for the temperature distribution in the melt, which are in a good agreement with the experimental data. This shows the opportunity of using this approach for engineering requirements. Numerical studies with different grid and solver parameters based on the unstationary 3D model show that numerical tools based on k- e turbulence model do not simulate the long period oscillations of the averaged melt flow. First results using LES turbulence model show low-frequency instabilities of the melt as carried out experimentally.


  1. E. Baake, A. M ühlbauer, B. Nacke: Heat and mass transfer in the turbulent melt flow of the crucible inductor furnace. International Colloquium on Modeling of Material Processing, (Riga, Latvia, 1999), pp.98-103
  2. A. Mühlbauer, E. Baake, A. Jakovics: Influence of the filling level on the turbulent material exchange in the melt of induction crucible furnaces. Second International Conference on Energy Transfer in Magnetohydrodynamic Flows, (Aussois, France, 1994), Vol.1, pp.349-358
  3. L. Meinecke: Numerische Simulation der Schmelzstr ömung und Temperaturverteilung in Tiegelinduktor öfen. Diplomarbeit, Universität Hannover, Hannover, 1999
  4. E. Baake, A. M ühlbauer, A. Jakovitsch, W. Andree: Extension of the k- e model for the numerical simulation of the melt flow in induction crucible furnaces. Metallurgical and Material Transactions, Vol. 26B, (1994), pp.529-535
  5. E. Baake: Grenzleistungs- und Aufkohlungsverhalten von Induktions-Tiegelöfen. VDI-Verlag, Düsseldorf, 1994
  6. L.D. Landau, J.M. Lifschitz: Gidrodinamika, Moscow, 1986, pp.65-68
  7. B.E. Launder, D.B. Spalding: The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering 3 (1974), pp.269-289.
  8. A. Bojarevics, V. Bojarevics, J. Gelfgat, K. Pericleous: Liquid metal turbulent flow dynamics in a cylindrical container with free surface: experiment and numerical analysis. International Colloquium on Modeling of Material Processing, (Riga, Latvia, 1999), pp.68-73