Numerical Simulation of the Flow and Heat Transfer in an Electric Steel Tempering Furnace

: Heat treatments, such as steel tempering, are temperature-controlled processes. It allows ferrous steel to stabilize its structure after the heat treatment and quenching stages. The tempering temperature also determines the hardness of the steel, preferably to its optimum working strength. In a tempering furnace, a heat-resistant fan is commonly employed to generate moderate gas circulation to obtain adequate temperature homogeneity and heat transfer. Nevertheless, there is a tradeo ﬀ because the overall thermal e ﬃ ciency is expected to reduce because of the high rotating speed of the fan. Therefore, this study numerically investigates the thermal e ﬃ ciency changes of an electric tempering furnace due to changes in the rotating speed of the fan and the e ﬀ ects on temperature homogeneity and the heat transfer rate to the load. Heat losses through the walls were calculated from the external temperature measurement of the furnace. Four di ﬀ erent speeds were simulated: 720, 990, 1350, and 1800 rpm. Thermal homogeneity was improved at higher rotating speeds; this is because the recirculation zone caused by the fan improved the ﬂow mixing and the heat transfer. However, it was found that the thermal e ﬃ ciency of the tempering furnace decreased as the rotating speed values increased. Therefore, these characteristics should be modulated to obtain a proﬁt when controlling the rotating speed. For example, although thermal e ﬃ ciency decreases by 20% when the rotating speed is doubled, the heat transfer rate to load is increased by up to 50%, which can be beneﬁcial in decreasing the process of tempering times. ﬁxed as the setpoint for the heating chamber. One of the aims of this stage is to remove any residual humidity on the steel parts. After two hours of preheating, the cooling stage starts. The setpoint is ﬁxed at 773 K, which is a temperature where the tempering stress relief starts [5]. After one hour, the setpoint is gradually increased to 20 K for half an hour. This process is repeated two times to obtain adequate steel hardness.


Introduction
Industrial furnaces have become vitally important equipment since they are involved in the production of many consumer products, such as food, beverages, containers, machining tools, infrastructure materials, amongst other applications [1]. In such furnaces, heat generation technology is mainly based on gas heating as well as electricity. Operation costs related to energy consumption are significant in most high-temperature furnaces; therefore, there is potential for energy savings [2,3].
Furnaces should uniformly perform their heat transfer process. Non-uniform temperature distribution inside the heating chamber may lead to low-quality products, thereby increasing production costs. In some applications, such as the heat treatment of metals, it is of utmost importance that the temperature is uniform inside the furnace [4,5].
For example, in a typical hardening process, steel parts are heated to a preset temperature between from 1000 K to 1470 K. Then, steel is quenched in special dielectric oil or by other means, such as molten salts or polymers, to rapidly remove heat. Quench-hardened steel is usually too brittle and must be tempered. The main goal of tempering is to decrease hardness and to increase the toughness of the ferrous material by relieving the tension that accumulates in the lattice of the metal [4,5]. This is unloading of parts, which is done by lifting the lid away and hooking the parts out of the furnace. The electric furnace is composed of two chambers: the heating chamber, where 22 metal coils are electrically heated, and the process chamber, where the steel parts are stored and tempered. The heating chamber is internally covered with insulating refractory bricks, classified as group 23 according to the ASTM C155 standard [14]. The process chamber is a cylinder with a 0.3 m radius and 1.2 m height. This space is covered with refractory concrete and AISI 304 stainless steel plates with a thickness of about 7.94 mm (5/16 inches) for mechanical and thermal protection. The gross load capacity is rated at 400 kg. Such a load can be put into metallic baskets and inserted from the top of the process chamber. and 1.2 m height. This space is covered with refractory concrete and AISI 304 stainless steel plates with a thickness of about 7.94 mm (5/16 inches) for mechanical and thermal protection. The gross load capacity is rated at 400 kg. Such a load can be put into metallic baskets and inserted from the top of the process chamber.
The furnace is powered by a 220 V three-phase system at a 60 Hz frequency, with a total nominal power of up to 50 kW and a maximum operating temperature of 1120 K, while the temperature in the heating coil can exceed 1370 K. The furnace has a fan at the bottom of the heating chamber that operates at a fixed rotating speed of 990 rpm, powered by a 2.24 kW (3 hp) electric motor. The fan was manufactured using AISI 310 stainless steel. The fan moves hot gases from the hot coils downward inside the heating chamber and then forces them upwards into the workspace and through the load. Then, gases are recirculated through the coils bank. The recirculated gas inside the equipment is a fixed mass of air, with no inlets or outlets. This furnace has two thermocouples to control the temperature: one located at the entrance of the gases to the heating chamber, near the top of the furnace (controlling the process temperature). The second thermocouple is located at the heating chamber itself, and its main function is to avoid coil overheating. Both inputs work together in controlling the heating power delivered to the metallic load.  The furnace is powered by a 220 V three-phase system at a 60 Hz frequency, with a total nominal power of up to 50 kW and a maximum operating temperature of 1120 K, while the temperature in the heating coil can exceed 1370 K. The furnace has a fan at the bottom of the heating chamber that operates at a fixed rotating speed of 990 rpm, powered by a 2.24 kW (3 hp) electric motor. The fan was manufactured using AISI 310 stainless steel. The fan moves hot gases from the hot coils downward inside the heating chamber and then forces them upwards into the workspace and through the load. Then, gases are recirculated through the coils bank. The recirculated gas inside the equipment is a fixed mass of air, with no inlets or outlets. This furnace has two thermocouples to control the temperature: one located at the entrance of the gases to the heating chamber, near the top of the furnace (controlling the process temperature). The second thermocouple is located at the heating chamber itself, and its main function is to avoid coil overheating. Both inputs work together in controlling the heating power delivered to the metallic load. Figure 2 presents the temporal evolution of the temperature setpoint inside the furnace during the tempering process studied in this work. This heat treatment follows a three-stage process: preheating, cooling, and sustaining. In the preheating stage, the furnace lid is opened, and a loaded charging basket is introduced using an overhead crane. Afterward, the furnace is closed, and a temperature of 973 K is Energies 2020, 13, 3655 4 of 22 fixed as the setpoint for the heating chamber. One of the aims of this stage is to remove any residual humidity on the steel parts. After two hours of preheating, the cooling stage starts. The setpoint is fixed at 773 K, which is a temperature where the tempering stress relief starts [5]. After one hour, the setpoint is gradually increased to 20 K for half an hour. This process is repeated two times to obtain adequate steel hardness.
in controlling the heating power delivered to the metallic load. Figure 2 presents the temporal evolution of the temperature setpoint inside the furnace during the tempering process studied in this work. This heat treatment follows a three-stage process: preheating, cooling, and sustaining. In the preheating stage, the furnace lid is opened, and a loaded charging basket is introduced using an overhead crane. Afterward, the furnace is closed, and a temperature of 973 K is fixed as the setpoint for the heating chamber. One of the aims of this stage is to remove any residual humidity on the steel parts. After two hours of preheating, the cooling stage starts. The setpoint is fixed at 773 K, which is a temperature where the tempering stress relief starts [5]. After one hour, the setpoint is gradually increased to 20 K for half an hour. This process is repeated two times to obtain adequate steel hardness.

Numerical Models
Steady-state operation of the furnace at the end of the tempering process is considered. It corresponds to a turbulent flow occurring in a low Mach number regime. The transport equations governing fluid flow for this kind of problem are the Favre-averaged equations transport equations [15]. First, the continuity equation is where ̅ and ̃ are the average density and velocity. The momentum equation is also solved: where , ̅ and g are laminar viscosity, averaged pressure, and gravity, respectively. In this equation, is the identity matrix. The Reynolds stresses −∇ • ( ̅ ′ ′ ) is a term that must be modeled. In this work, the Boussinesq hypothesis was followed:

Numerical Models
Steady-state operation of the furnace at the end of the tempering process is considered. It corresponds to a turbulent flow occurring in a low Mach number regime. The transport equations governing fluid flow for this kind of problem are the Favre-averaged equations transport equations [15]. First, the continuity equation is ∂ρ ∂t where ρ and v are the average density and velocity. The momentum equation is also solved: where µ l , p and g are laminar viscosity, averaged pressure, and gravity, respectively. In this equation, I is the identity matrix. The Reynolds stresses −∇· ρ v v is a term that must be modeled. In this work, the Boussinesq hypothesis was followed: where k is the averaged turbulent kinetic energy: The turbulent viscosity µ t is obtained from the following equation: Energies 2020, 13, 3655 where ε is the averaged dissipation of turbulent kinetic energy and C µ is a modeling constant equal to 0.09, taken from the standard k-epsilon turbulence model. This turbulence model has been used successfully in the past for simulations of thermal treatment furnaces [16]. In this model, two additional transport equations for k and ε must be solved: where the model constants C ε1 , C ε2 , σ k , σ ε are 1.44, 1.92, 1, and 1.3, respectively. Additionally, standard wall functions for modeling the velocity profile near the walls were employed in this work [17]. The energy transport equation is also solved in this work: where k t is the thermal conductivity of the air, T is the average temperature and S rad is the source term due to radiation. The total non-chemical energy E is defined as where h in this case stands for enthalpy. In this work, radiation was modeled using a discrete ordinate (DO) model. The DO model solves the radiative heat transfer equation (RTE) for a finite number of discrete angles, each one associated with a vector direction → s in a polar system (θ, ϕ, z): where I → s is the spectral radiation intensity, a is the absorption coefficient, σ s is the dispersion coefficient, Ω is the solid angle, φ is the scattering phase function and σ is the Stefan-Boltzmann constant (5.67 × 10 −8 Wm −2 K −4 ). This radiation model requires θ and ϕ angles to be discretized in a finite number of subdivisions; the more subdivisions, the better the results, but the computational time increases substantially. In this work, 2 × 2 subdivisions were used. Additionally, because isotropic scattering was assumed in this work, the scattering phase function φ → s , → s in the RTE is unity. Finally, the source term in the energy equation is [18]: where n is the refractive index and G is the incident radiation, which is evaluated as IdΩ.
The baseline operation simulation considers the rotating speed of the fan to be 990 rpm. Simulations at 720, 1350, and 1800 rpm were also performed. The fan is configured as a rotating wall using a moving-reference frame (MRF) method with a frozen rotor approach-as presented in Figure 3. MRF is an acceptable modeling approach that has been used successfully in axial [19], centrifugal [20], and other kinds of fans [21]. Additionally, a previous study using a nitrocarburation furnace with a fan [22] was modeled using the MRF method.

Boundary Conditions
To properly configure boundary conditions in the simulations, the heat flux lost through each wall was estimated with the experimentally measured external wall temperature. Temperatures were measured using a type K thermocouple of 1 mm point diameter with a UT320 UNI-T digital indicator. Measurements were made at 24 points on the outside of the furnace, as shown in Figure 4. This figure specifies the location of each temperature measuring point using a dot. The temperature was measured three times at each point and then they were averaged for the whole surface. The maximum standard deviation of the measurement for any point was 12.4 K, while the average standard deviation was 4.5 K.

Boundary Conditions
To properly configure boundary conditions in the simulations, the heat flux lost through each wall was estimated with the experimentally measured external wall temperature. Temperatures were measured using a type K thermocouple of 1 mm point diameter with a UT320 UNI-T digital indicator. Measurements were made at 24 points on the outside of the furnace, as shown in Figure 4. This figure specifies the location of each temperature measuring point using a dot. The temperature was measured three times at each point and then they were averaged for the whole surface. The maximum standard deviation of the measurement for any point was 12.4 K, while the average standard deviation was 4.5 K.
Boundary conditions are presented in Table 1. The domain simulated in this work corresponds to a closed system with no inlets or outlets. Due to its symmetry, only half of the fluid domain was simulated. The heat flux presented in Table 1 was calculated from temperature measurements and was then configured as the heat transfer boundary condition for each outer wall of the furnace.  q was assumed to be due to a combination of convection and radiation: .
The tempering furnace studied in this work is located inside a room with no significant airflow coming in or out. Therefore, a 5 Wm −2 K −1 natural convective coefficient h conv was assumed [23], and the convective heat flux . q conv was calculated with Newton's law of cooling: Energies 2020, 13, 3655 where the surrounding temperature T surr is 300 K on average and T s is the external wall surface temperature. Ts were calculated as the temperature average of every section of the oven ( Figure 4): side heating chamber (4 data), back heating chamber (2 data), workspace chamber (6 data) and, top heating chamber (12 data). To obtain a measurement of the radiation losses from the walls, the heat flux was calculated with the following expression [24]: .
where ε is the surface emissivity. All surfaces were considered opaque with a diffusive fraction equal to 1. Refractory emissivity was assumed to be 0.65, and steel emissivity was assumed to be 0.85. Air was modeled as an ideal gas with variable properties depending on temperatures, such as density or viscosity. Pressure-momentum coupling was achieved using a coupled algorithm for faster convergence. Convective terms in all transport equations had second-order accuracy and convergence criteria were established at 1 × 10 −6 for all residuals. Simulations were carried out in ANSYS FLUENT ® (Canonsburg, PA, USA); v19 in a computer equipped with an Intel Core i7 8700k ® processor with 16 GB of RAM.  Boundary conditions are presented in Table 1. The domain simulated in this work corresponds to a closed system with no inlets or outlets. Due to its symmetry, only half of the fluid domain was simulated. The heat flux presented in Table 1 was calculated from temperature measurements and was then configured as the heat transfer boundary condition for each outer wall of the furnace.

Mesh Independence Study
The geometry was constructed in a the-dimensional computational domain, which comprised components such as a heating chamber, fan, metal coils, and process chamber. The geometry was taken from a real pit type furnace. The metal coils were constructed as cylindrical metal bars. Additionally, the present study constructed non-structured full tetrahedral meshes using ANSYS Meshing ® (Canonsburg, PA, USA), as shown in Figure 5. After a grid-independent test, the coarse mesh with 1.5 million cells was rejected, due to the simulations not capturing the thermal behavior of the furnaces well, mainly near the wall zones.
ANSYS FLUENT ® (Canonsburg, PA, USA); v19 in a computer equipped with an Intel Core i7 8700k ® processor with 16 GB of RAM.

Mesh Independence Study
The geometry was constructed in a three-dimensional computational domain, which comprised components such as a heating chamber, fan, metal coils, and process chamber. The geometry was taken from a real pit type furnace. The metal coils were constructed as cylindrical metal bars. Additionally, the present study constructed non-structured full tetrahedral meshes using ANSYS Meshing ® (Canonsburg, PA, USA), as shown in Figure 5. After a grid-independent test, the coarse mesh with 1.5 million cells was rejected, due to the simulations not capturing the thermal behavior of the furnaces well, mainly near the wall zones.  Two other meshes were considered in this work: one fine-mesh approximately double the size of the coarse mesh (3 million cells) and a very fine mesh composed of about 6 million cells. Figure 6 shows the results of the temperature as a function of the radial distance of the workspace or process chamber of the furnace (at Z = 0.5 m). This horizontal line is located at the middle of the vertical symmetry plane passing through the center of the process chamber. In this case, the baseline operation for the three meshes shows that there is less than a 2% difference in the results between the fine and the very fine meshes. Therefore, it can be argued that the fine mesh is accurate enough for this study and the 3 million-cell mesh is chosen as the mesh for all the results that are reported in the next section. This mesh represents a good compromise between detail and computational time. Simulations are initialized at 700 K and convergence is reached at close to 10 thousand iterations. The computational time for each run was close to 100 h. Two other meshes were considered in this work: one fine-mesh approximately double the size of the coarse mesh (3 million cells) and a very fine mesh composed of about 6 million cells. Figure 6 shows the results of the temperature as a function of the radial distance of the workspace or process chamber of the furnace (at Z = 0.5 m). This horizontal line is located at the middle of the vertical symmetry plane passing through the center of the process chamber. In this case, the baseline operation for the three meshes shows that there is less than a 2% difference in the results between the fine and the very fine meshes. Therefore, it can be argued that the fine mesh is accurate enough for this study and the 3 million-cell mesh is chosen as the mesh for all the results that are reported in the next section. This mesh represents a good compromise between detail and computational time.
Simulations are initialized at 700 K and convergence is reached at close to 10 thousand iterations. The computational time for each run was close to 100 h. Figure 6. Temperature profile at a horizontal line located in the middle of the workspace for the three different meshes studied at 990 rpm.

Convective Coefficient to the Load and Energy Balance
Simulations were performed with and without the metallic charging basket needed for hanging or supporting the steel parts to be treated inside the process chamber of the furnace. An empty and fully loaded basket are presented in Figure 7. Such a basket could affect the volumetric flow supplied

Convective Coefficient to the Load and Energy Balance
Simulations were performed with and without the metallic charging basket needed for hanging or supporting the steel parts to be treated inside the process chamber of the furnace. An empty and fully loaded basket are presented in Figure 7. Such a basket could affect the volumetric flow supplied by the fan, and therefore, its influence on air-flow was analyzed as well. The basket is made of AISI 304 stainless steel, with thermal conductivity of 16.6 W/(m K), a density of 7900 kg/m 3 , and specific heat of 477 J/(kg K) [23]. Calculated from density and volume, the mass of this basket is 54.7 kg. In this work, the heating of a cylindrical steel pin, which is a typical part of this kind of furnace and process, was considered. Each pin is 300 mm long with a 50.8 mm external and 25.4 mm internal diameter, made of carbon steel, with thermal conductivity of 56.7 W/(m K), a density of 7854 kg/m 3, and specific heat of 487 J/(kg K) [23]. Overall, the heating of 180 kg of these pins was considered in this work.
The results predicted by simulations allowed for the calculation of heat transfer to the load by convection using three well-known empirical correlations for the average non-dimensional Nusselt number Nu. The maximum value of axial speed predicted by the CFD software, as well as the process temperature, were the main inputs for these correlations. Then, the average of the correlations is reported in the results. The first correlation considered in this work is the one developed by Hilpert [23,25]: where parameters C and m are obtained from Reference [23] and they depend on the Reynolds number Re, with the pin diameter as the characteristic length. Prandtl number Pr is obtained with the surface-fluid average temperature, using air as a fluid. The second correlation used was developed by Žukauskas [26]: where parameters C and m are obtained from Reference [23] and the Prs is the Prandtl number based on the surface temperature, considered in this work to be 773.15 K. In addition to these correlations, the expression of Churchill and Bernstein [27], where all properties are evaluated at the temperature of the film, was also considered: In this work, the heating of a cylindrical steel pin, which is a typical part of this kind of furnace and process, was considered. Each pin is 300 mm long with a 50.8 mm external and 25.4 mm internal diameter, made of carbon steel, with thermal conductivity of 56.7 W/(m K), a density of 7854 kg/m 3 , and specific heat of 487 J/(kg K) [23]. Overall, the heating of 180 kg of these pins was considered in this work.
The results predicted by simulations allowed for the calculation of heat transfer to the load by convection using three well-known empirical correlations for the average non-dimensional Nusselt number Nu. The maximum value of axial speed predicted by the CFD software, as well as the process temperature, were the main inputs for these correlations. Then, the average of the correlations is reported in the results. The first correlation considered in this work is the one developed by Hilpert [23,25]: where parameters C and m are obtained from Reference [23] and they depend on the Reynolds number Re, with the pin diameter as the characteristic length. Prandtl number Pr is obtained with the surface-fluid average temperature, using air as a fluid. The second correlation used was developed by Žukauskas [26]: Energies 2020, 13, 3655 10 of 22 where parameters C and m are obtained from Reference [23] and the Pr s is the Prandtl number based on the surface temperature, considered in this work to be 773.15 K. In addition to these correlations, the expression of Churchill and Bernstein [27], where all properties are evaluated at the temperature of the film, was also considered: The energy balance for a closed furnace during the heating stage has recently been studied in Reference [28]. In this case, energy input is not due to the combustion of natural gas, but rather electricity, and thermal efficiency η can therefore be calculated as Heat trans f er rate to the load Electric power input , where the heat transfer to the load corresponds to the sensible heat demanded by the basket and the cylindrical steel pins during a processing time of 2 h, and the electric power input is calculated with the energy required by the fan and the metal coils. In theory, the fan studied in this work should follow the similarity law [24]: where .
Q. is the volumetric flow rate, ω is the angular velocity, and D is the fan's outer diameter. This equation states that, for any two different operational conditions 1 and 2, the dimensionless flow coefficient Φ is assumed to be equal. If the fans are also assumed to have the same dimensionless power coefficient C P , then the following similarity law also should apply: where . W sha f t is the shaft power. These laws are derived from dimensional analysis following Buckingham's theorem π. An analysis of the validity of these similarity laws of the fan modeled in this work is presented in the results.

Results and Discussion
The numerical results include a baseline operation and the effect of the fan rotating speed in the tempering furnace. We focus on the temperature and fluid dynamic fields inside the furnace, as well as the heat transfer behavior and efficiencies when the speed is changed.

Baseline Operation
The tempering system is characterized by its capacity to achieve a homogeneous temperature field throughout the process chamber, with the help of an internal fan. Figure 8 shows the temperature distribution at the vertical symmetry plane passing through the center of the heating and the process chambers of the tempering furnace. This corresponds to a baseline operation that has a fan rotating speed of 990 rpm. From the figure, it can be seen that the global thermal behavior of both chambers is well captured by the numerical calculations: the highest temperatures are obtained at the metal coils, the heat transfer is forced by the fan from the higher temperature zone to the lower temperature zone, and the thermal homogeneity occurs in the process chamber. Additionally, the experimental temperature data inside the furnace is in agreement with the numerical data. For example, the temperature of the process chamber was 813.1 K (last stage of Figure 2) and the simulated temperatures correspond to 806.4 K.
Energies 2020, 13, x FOR PEER REVIEW 12 of 23 through the walls, and therefore was not taken into account (the adiabatic wall was configured as presented in Table 1).
(a) (b)  Figure 9 presents streamlines for the baseline operation without a charging basket. Flow velocity near the rotor reaches up to 16.8 m/s. In contrast, velocity inside the process chamber is considerably lower, less than 2 m/s in most of the chamber. Despite the furnace being of convective type, this result indicates that the actual velocity, in the zone where the steel parts are exposed to hot air, is relatively low. As was presented in the methodology, convection heat transfer is proportional to the convective coefficient. This coefficient increases with the non-dimensional Reynolds number, which in turn is proportional to the fluid velocity. A low velocity within the workspace means a low heat transfer rate to the load and higher processing times.
Streamlines in Figure 9 also show an important recirculation zone inside the process chamber of the furnace. Recirculation helps increase the residence time of the hot gas in the workspace. This inner swirl has a similar function as the inner baffles reported by other authors [1] because it helps to augment the time that the hot gases remain in the process chamber of the furnace. In addition, some eddies are shown at the upper corners of the heating and process chambers, which could hinder the fluid transport inside the furnace. Unlike the main recirculation zone, these eddies occur in regions where they are not relevant to the process.
Besides, from Figure 9 it can be seen that a significant amount of fluid coming from the process chamber is directed to the left side wall of the metal coils, instead of the metal coils. This accords with the original design of the furnace, which does not allow for optimally guiding the fluid to the coils. A geometrical modification at the upper heating chamber zone could improve the heat transfer from the metal coils to the recirculating air. The temperatures obtained inside the workspace chamber are between 750 K and 800 K. Therefore, the temperature uniformity range of the furnace is about 50 K, according to the DIN 17052-1 standard. Nevertheless, the four cut views presented at different heights in Figure 8 show that, as fluid flows from the bottom to the top of the chamber, the outer temperature next to the walls becomes closer to the center level and can be regarded as having a better temperature distribution, albeit equal in uniformity according to the norm. The results also show that the rotor is exposed to about 800 K during operation. The shaft of this rotor must be actively cooled to avoid heat damage to the engine; however, this heat loss was considered to be negligible, compared to the heat losses through the walls, and therefore was not taken into account (the adiabatic wall was configured as presented in Table 1). Figure 9 presents streamlines for the baseline operation without a charging basket. Flow velocity near the rotor reaches up to 16.8 m/s. In contrast, velocity inside the process chamber is considerably lower, less than 2 m/s in most of the chamber. Despite the furnace being of convective type, this result indicates that the actual velocity, in the zone where the steel parts are exposed to hot air, is relatively low. As was presented in the methodology, convection heat transfer is proportional to the convective coefficient. This coefficient increases with the non-dimensional Reynolds number, which in turn is proportional to the fluid velocity. A low velocity within the workspace means a low heat transfer rate to the load and higher processing times.
Streamlines in Figure 9 also show an important recirculation zone inside the process chamber of the furnace. Recirculation helps increase the residence time of the hot gas in the workspace. This inner swirl has a similar function as the inner baffles reported by other authors [1] because it helps to augment the time that the hot gases remain in the process chamber of the furnace. In addition, some eddies are shown at the upper corners of the heating and process chambers, which could hinder the fluid transport inside the furnace. Unlike the main recirculation zone, these eddies occur in regions where they are not relevant to the process.
Besides, from Figure 9 it can be seen that a significant amount of fluid coming from the process chamber is directed to the left side wall of the metal coils, instead of the metal coils. This accords with the original design of the furnace, which does not allow for optimally guiding the fluid to the coils. A geometrical modification at the upper heating chamber zone could improve the heat transfer from the metal coils to the recirculating air. The occurrence of a wide recirculation zone is evident in Figure 9. The figure also shows the mass recirculation at several heights of the process chamber (Z-direction). This was calculated taking into account the principle of mass conservation. That is, the mass flux in the Z-direction of the process chamber remains constant. This fact allows the use of negative air velocities to calculate the air mass recirculation in different transversal planes and establish a profile of recycled gases that increases progressively until reaching a maximum near the center of the vortex and then diminishes (right side of Figure 9). This result is similar to the one found by a previous study [29] using a non-reactive coflow jet. This type of calculation has been used with reactive gases [30,31], providing valuable fluid dynamics information of combustion chambers. In addition, it can be observed that even at high Zdirection of the process chamber, there is an air mass recirculation. This can be explained due to the airflow hits the top wall of the camera, and then come back.
Prediction of the behavior of the fluid near the walls is relevant, as shown above. Figure 10 displays the distribution of the Y+ on the inner walls of the furnace. Y+ is a non-dimensional measure of the first element size on the wall of the computational mesh. The maximum value of Y+ is about 65 on the walls next to the fan, which is lower than 100 and can be considered an acceptable value for the standard wall function of the k-epsilon turbulence model [17]. Therefore, this result confirms that the mesh resolution is adequate near the walls. The occurrence of a wide recirculation zone is evident in Figure 9. The figure also shows the mass recirculation at several heights of the process chamber (Z-direction). This was calculated taking into account the principle of mass conservation. That is, the mass flux in the Z-direction of the process chamber remains constant. This fact allows the use of negative air velocities to calculate the air mass recirculation in different transversal planes and establish a profile of recycled gases that increases progressively until reaching a maximum near the center of the vortex and then diminishes (right side of Figure 9). This result is similar to the one found by a previous study [29] using a non-reactive coflow jet. This type of calculation has been used with reactive gases [30,31], providing valuable fluid dynamics information of combustion chambers. In addition, it can be observed that even at high Z-direction of the process chamber, there is an air mass recirculation. This can be explained due to the airflow hits the top wall of the camera, and then come back.
Prediction of the behavior of the fluid near the walls is relevant, as shown above. Figure 10 displays the distribution of the Y+ on the inner walls of the furnace. Y+ is a non-dimensional measure of the first element size on the wall of the computational mesh. The maximum value of Y+ is about 65 on the walls next to the fan, which is lower than 100 and can be considered an acceptable value for the standard wall function of the k-epsilon turbulence model [17]. Therefore, this result confirms that the mesh resolution is adequate near the walls. Figure 11 shows the temperature distribution in the symmetry boundary of the furnace for the baseline operation (990 rpm) with the metallic charging basket inside. In this case, the temperature uniformity range in the process chamber (the right side of the figure) is higher, at about 70 K; meanwhile, in the case without the charging basket, it is about 50 K (Figure 8). This homogeneity is quantified by applying the concept of temperature uniformity, which represents the difference in temperature between different points in the system. In this case, an acceptable calculation of temperature uniformity is obtained as the maximum value of the temperature difference between the central point of the furnace and the corner points [10]. As can be seen in the Figure, the gases near the left lateral wall of the process chamber are cooler than the rest of the gases in the chamber. Besides, the upper part of the basket is cooler than the lower part. In a real process, treatments with noticeable differences in temperature within the chamber could impact the optimum steel working strength. Therefore, strategies to improve the thermal homogeneity of the furnace (e.g., fan rotating speed variation, as proposed in this study) should be developed.  Figure 11 shows the temperature distribution in the symmetry boundary of the furnace for the baseline operation (990 rpm) with the metallic charging basket inside. In this case, the temperature uniformity range in the process chamber (the right side of the figure) is higher, at about 70 K; meanwhile, in the case without the charging basket, it is about 50 K (Figure 8). This homogeneity is quantified by applying the concept of temperature uniformity, which represents the difference in temperature between different points in the system. In this case, an acceptable calculation of temperature uniformity is obtained as the maximum value of the temperature difference between the central point of the furnace and the corner points [10]. As can be seen in the Figure, the gases near the left lateral wall of the process chamber are cooler than the rest of the gases in the chamber. Besides, the upper part of the basket is cooler than the lower part. In a real process, treatments with noticeable differences in temperature within the chamber could impact the optimum steel working strength. Therefore, strategies to improve the thermal homogeneity of the furnace (e.g., fan rotating speed variation, as proposed in this study) should be developed.   Figure 11 shows the temperature distribution in the symmetry boundary of the furnace for the baseline operation (990 rpm) with the metallic charging basket inside. In this case, the temperature uniformity range in the process chamber (the right side of the figure) is higher, at about 70 K; meanwhile, in the case without the charging basket, it is about 50 K (Figure 8). This homogeneity is quantified by applying the concept of temperature uniformity, which represents the difference in temperature between different points in the system. In this case, an acceptable calculation of temperature uniformity is obtained as the maximum value of the temperature difference between the central point of the furnace and the corner points [10]. As can be seen in the Figure, the gases near the left lateral wall of the process chamber are cooler than the rest of the gases in the chamber. Besides, the upper part of the basket is cooler than the lower part. In a real process, treatments with noticeable differences in temperature within the chamber could impact the optimum steel working strength. Therefore, strategies to improve the thermal homogeneity of the furnace (e.g., fan rotating speed variation, as proposed in this study) should be developed.   Figure 12 presents streamlines for the baseline operation with a charging basket. The maximum velocity is about 17 m/s in the rotor zone, while velocities inside and around the basket are closer to 2 m/s in most of the chamber. It should be noted that the inner recirculation zone identified in the case without a basket is still present. Additionally, a significant amount of fluid coming from the process chamber, directed to the left side wall of the metal coils, is obtained. Therefore, velocity values, as well as distributions, are very similar in cases with and without a metallic charging basket. This is due to the basket geometry configuration allowing for the flow of gases to pass through it and it hardly offering axial resistance. An analogy can be drawn with the temperature case profile, where some differences in thermal homogeneity were found. In that case, the basket offers thermal resistance during the preheating stage, and there are energy losses through the walls during the steady-state stage. This is due to the basket geometry configuration allowing for the flow of gases to pass through it and it hardly offering axial resistance. An analogy can be drawn with the temperature case profile, where some differences in thermal homogeneity were found. In that case, the basket offers thermal resistance during the preheating stage, and there are energy losses through the walls during the steady-state stage.  Figure 13 presents energy inflows and outflows for the furnace for the baseline operation of 990 rpm, using two Sankey diagrams-one during the preheating stage and the other at the steady-state. It should be noted that input energy demand is at its maximum during the preheating stage. When the load reaches the steady-state condition, electricity input diminishes to the minimum level required to maintain the temperature inside the process chamber, because heat is being lost through the walls. In this furnace, the energy loss through the walls is roughly half of the energy required for heating the load, suggesting that furnace insulation should be improved. It should also be noted that heat recirculation of hot air between chambers accounts for up to 83.7% of all energy involved in the system. This flow is generated by the fan at a low energy expense of 1.34 kW.
If the rotating speed is increased from 990 rpm, fan energy demand would also increase, and efficiency would, therefore, be decreased. However, the heat transfer rate to the load could be improved, and shorter tempering times could be achieved. As pointed out in Reference [5], steel quality is not negatively affected by longer tempering times, but shorter processing times are desirable in a furnace that operates continuously. Therefore, the heat transfer rate could benefit from increased rotating speed. This behavior is studied and presented in the following section.  Figure 13 presents energy inflows and outflows for the furnace for the baseline operation of 990 rpm, using two Sankey diagrams-one during the preheating stage and the other at the steady-state. It should be noted that input energy demand is at its maximum during the preheating stage. When the load reaches the steady-state condition, electricity input diminishes to the minimum level required to maintain the temperature inside the process chamber, because heat is being lost through the walls. In this furnace, the energy loss through the walls is roughly half of the energy required for heating the load, suggesting that furnace insulation should be improved. It should also be noted that heat recirculation of hot air between chambers accounts for up to 83.7% of all energy involved in the system. This flow is generated by the fan at a low energy expense of 1.34 kW.

Effect of the Fan Rotating Speed
To estimate the temperature uniformity induced by the flow mixing in the furnace, the temperature distributions were analyzed. The temperature distribution for different fan rotating speeds at a horizontal line located in the middle of the process chamber (at Z = 0.5 m) is presented in If the rotating speed is increased from 990 rpm, fan energy demand would also increase, and efficiency would, therefore, be decreased. However, the heat transfer rate to the load could be improved, and shorter tempering times could be achieved. As pointed out in Reference [5], steel quality is not negatively affected by longer tempering times, but shorter processing times are desirable in a furnace that operates continuously. Therefore, the heat transfer rate could benefit from increased rotating speed. This behavior is studied and presented in the following section.

Effect of the Fan Rotating Speed
To estimate the temperature uniformity induced by the flow mixing in the furnace, the temperature distributions were analyzed. The temperature distribution for different fan rotating speeds at a horizontal line located in the middle of the process chamber (at Z = 0.5 m) is presented in Figure 14 for the case of an empty furnace. The higher the rotating speed, the higher the average temperature inside the chamber. This trend could be explained because the tempering furnace was analyzed as a closed system with no inlets or outlets. As the rotating speed increases, the required fan power also increases; therefore, as the losses through the walls were similar to each other, the internal energy inside the chamber, and the temperature increased.

Effect of the Fan Rotating Speed
To estimate the temperature uniformity induced by the flow mixing in the furnace, the temperature distributions were analyzed. The temperature distribution for different fan rotating speeds at a horizontal line located in the middle of the process chamber (at Z = 0.5 m) is presented in Figure 14 for the case of an empty furnace. The higher the rotating speed, the higher the average temperature inside the chamber. This trend could be explained because the tempering furnace was analyzed as a closed system with no inlets or outlets. As the rotating speed increases, the required fan power also increases; therefore, as the losses through the walls were similar to each other, the internal energy inside the chamber, and the temperature increased.
From Figure 14, the maximum temperature difference of the baseline operation of 990 rpm, between radial distances of -0.25 and 0.25 m, is 47.5 K. This difference is reduced to 41.2 and 34.6 K at 1350 and 1800 rpm, respectively. Therefore, a higher rotating speed of the fan promotes temperature homogeneity. Additionally, it should be noted that the temperature profile is not symmetric around the process chamber centerline, and in all cases, the temperature decreases when the radial distance is very close to the walls (0.3 or -0.3 m). Figure 14. Temperature profile at a horizontal line located in the middle of the workspace for the five different angular rotation rates studied. Figure 15 shows the axial velocity profile at a horizontal line located in the middle of the workspace for the four different rotating speeds studied (same line presented in Figure 14). Flow is From Figure 14, the maximum temperature difference of the baseline operation of 990 rpm, between radial distances of −0.25 and 0.25 m, is 47.5 K. This difference is reduced to 41.2 and 34.6 K at 1350 and 1800 rpm, respectively. Therefore, a higher rotating speed of the fan promotes temperature homogeneity. Additionally, it should be noted that the temperature profile is not symmetric around the process chamber centerline, and in all cases, the temperature decreases when the radial distance is very close to the walls (0.3 or −0.3 m). Figure 15 shows the axial velocity profile at a horizontal line located in the middle of the workspace for the four different rotating speeds studied (same line presented in Figure 14). Flow is not symmetric around the chamber centerline and both positive and negative velocities appear. A negative velocity axial means that the flow is reversing from top to bottom, which is a result that, combined with the positive velocity values, indicates the presence of a swirl. Therefore, there is an inner recirculation zone within the process chamber due to a swirl that intensifies with higher rpms, and as axial velocities differentials become larger. As shown before, for the 1800 rpm case, temperature homogeneity is enhanced, and the velocity profiles suggest that this is due to a stronger swirl than in other cases. The fan thus plays an important role in enhancing the flow mixing and convective heat transfer in the furnace. The following analyses quantify the fan rotating speed effect in the convective heat transfer. and as axial velocities differentials become larger. As shown before, for the 1800 rpm case, temperature homogeneity is enhanced, and the velocity profiles suggest that this is due to a stronger swirl than in other cases. The fan thus plays an important role in enhancing the flow mixing and convective heat transfer in the furnace. The following analyses quantify the fan rotating speed effect in the convective heat transfer. The use of negative air velocities allowed the calculation of the mass recirculation (mr) within the process chamber. The results of this calculation can be seen in Figure 16a. A mass recirculation profile is displayed for each rotating speed which shows that as the gases proceed upwards (i.e., z increases) more gases are entrained into the upcoming gases until they reach a maximum located near the center of the vortex (see Figure 9). Two aspects are interesting to point out: as the rotating speed increases the mass recirculation also increases and the maximum moves to lower heights. In other words, more mass is recirculated, and the center of the vortex changes its Z-direction with the rotating speed increase. The maximum mr values are obtained in the range of 0.4 to 0.7 m. At the particular case of 720 rpm rotating speed, the maximum mr is located near a Z-direction of 0.7 m; and the mr values can be even higher than others rpms (i.e., 990 and 1350) in the upper zone of the process chamber. The use of negative air velocities allowed the calculation of the mass recirculation (mr) within the process chamber. The results of this calculation can be seen in Figure 16a. A mass recirculation profile is displayed for each rotating speed which shows that as the gases proceed upwards (i.e., z increases) more gases are entrained into the upcoming gases until they reach a maximum located near the center of the vortex (see Figure 9). Two aspects are interesting to point out: as the rotating speed increases the mass recirculation also increases and the maximum moves to lower heights. In other words, more mass is recirculated, and the center of the vortex changes its Z-direction with the rotating speed increase. The maximum mr values are obtained in the range of 0.4 to 0.7 m. At the particular case of 720 rpm rotating speed, the maximum mr is located near a Z-direction of 0.7 m; and the mr values can be even higher than others rpms (i.e., 990 and 1350) in the upper zone of the process chamber. The mass recirculation could be related to the thermal homogeneity inside the chamber. To do that we have used a temperature coefficient of uniformity (CU) following the definition presented in a previous study [16]. The coefficient was calculated in different transversal planes taking temperature data of each plane. This definition gives more information than the simple temperature difference presented in Section 3.1.
The CU profiles reported as a function of the Z-direction in Figure 16b evidence an increase in the temperature uniformity inside the process chamber as the rotating speed increases. It should be noted that even a low CU increase (i.e., 0.01) allows improvements in thermal homogeneity. For instance, in Figure 14 the maximum temperature difference at 990 rpm case was 47.5 K and the difference was reduced to 34.6 K in the 1800 rpm case. The temperature uniformity rise to high values indicates a sort of well-mixed condition, which is a valuable characteristic in convection type furnaces. In addition, the 720 rpm case is interesting because the vortex of the recirculation zone is located in the upper zone of the furnace. Therefore, it presents higher values of mr than other rotating The mass recirculation could be related to the thermal homogeneity inside the chamber. To do that we have used a temperature coefficient of uniformity (CU) following the definition presented in a previous study [16]. The coefficient was calculated in different transversal planes taking temperature data of each plane. This definition gives more information than the simple temperature difference presented in Section 3.1.
The CU profiles reported as a function of the Z-direction in Figure 16b evidence an increase in the temperature uniformity inside the process chamber as the rotating speed increases. It should be noted that even a low CU increase (i.e., 0.01) allows improvements in thermal homogeneity. For instance, in Figure 14 the maximum temperature difference at 990 rpm case was 47.5 K and the difference was reduced to 34.6 K in the 1800 rpm case. The temperature uniformity rise to high values indicates a sort of well-mixed condition, which is a valuable characteristic in convection type furnaces. In addition, the 720 rpm case is interesting because the vortex of the recirculation zone is located in the upper zone of the furnace. Therefore, it presents higher values of mr than other rotating speed conditions in that zone, according to Figure 16a. It also allows an improvement of temperature CU in that area, according to Figure 16b. This behavior suggests a correspondence between mass recirculation and thermal homogeneity.
Furthermore, some observations can be made from Figure 16b. At low Z-direction values, the fluid coming from the heating to the process chamber suffers a volume expansion that causes the temperature to drop; as a result, the temperature CU decreases. In Figures 9 and 14 can be seen that temperatures profiles are not symmetric around the process chamber centerline. This probably means that, even though the temperature values at the several locations become close to each other, some hot spots are present in some locations inside the chamber. This asymmetric temperature behavior decreases at higher heights of the chamber. The above is consistent with the rise of the temperature CU showed in Figure 16b.
The relative mass recirculation is presented in Figure 17. It normalizes the maximum mass recirculated (mr max ) over the total mass (mt). Maximum mr and mt increase with the rotating speed. However, the total mass rises at a higher rate. As a result, relative mass recirculation tends to decrease. For example, it takes a value of 0.87 at the 720 rpm case and 0.51 at 1800 rpm. These values mean that 87% or 51% of the inlet mass recirculates near to the vortex of the recirculation zone. These relative mass recirculation values are low when compared with others found in the literature [29][30][31] that use systems with high momentum coflow jets. The above suggests that the fluid dynamics characteristics of the tempering furnace could be improved by means of exploring other ways to increase the air momentum or even the use of secondary air. Energies 2020, 13, x FOR PEER REVIEW 19 of 23  Figure 18 presents the variation in the hot air flow rate inside the furnace with the fan rotating speed. Simulation results closely follow the first theoretical similarity law for fans, described by Equation (20), which assumes an equal dimensionless flow coefficient Φ for any two different operational conditions. This validates that the fan used in this study follows the classic laws for fans with low uncertainties.  Figure 18 presents the variation in the hot air flow rate inside the furnace with the fan rotating speed. Simulation results closely follow the first theoretical similarity law for fans, described by Equation (20), which assumes an equal dimensionless flow coefficient Φ for any two different operational conditions. This validates that the fan used in this study follows the classic laws for fans with low uncertainties. Figure 18 presents the variation in the hot air flow rate inside the furnace with the fan rotating speed. Simulation results closely follow the first theoretical similarity law for fans, described by Equation (20), which assumes an equal dimensionless flow coefficient Φ for any two different operational conditions. This validates that the fan used in this study follows the classic laws for fans with low uncertainties. Additionally, because the geometry of the fan remains constant, and the fluid is considered incompressible, the second similarity law, described in Equation (21), applies and assumes that the same power coefficient for any two different operational conditions. The power requirements are estimated in the next equation, which states that the power is proportional to the cube of the rotating speed:

Thermal Efficiency and Heat Transfer to the Load
Then, the calculated fan power takes the following values: 0.51, 1.34, 3.38, and 10.72 kW, with rotating speeds of 720, 990, 1350, and 1800 rpms, respectively. From the results, rotating speed variation can negatively affect the thermal efficiency of the furnace; however, it should be contrasted with its positives effect on thermal homogeneity and the convective heat transfer to the load. Figure 19 presents the heat transfer by convection to the load, represented by the dimensionless Nusselt number, and the thermal efficiency η as a function of the fan rotating speed. Both can be approximated to a quadratic equation as a function of the rpms with a coefficient of determination R 2 of 0.99. At the baseline operation of 990 rpm, the efficiency is 63.96%, while the Nusselt number is Additionally, because the geometry of the fan remains constant, and the fluid is considered incompressible, the second similarity law, described in Equation (21), applies and assumes that the same power coefficient C P for any two different operational conditions. The power requirements are estimated in the next equation, which states that the power is proportional to the cube of the rotating speed: Then, the calculated fan power takes the following values: 0.51, 1.34, 3.38, and 10.72 kW, with rotating speeds of 720, 990, 1350, and 1800 rpms, respectively. From the results, rotating speed variation can negatively affect the thermal efficiency of the furnace; however, it should be contrasted with its positives effect on thermal homogeneity and the convective heat transfer to the load. Figure 19 presents the heat transfer by convection to the load, represented by the dimensionless Nusselt number, and the thermal efficiency η as a function of the fan rotating speed. Both can be approximated to a quadratic equation as a function of the rpms with a coefficient of determination R 2 of 0.99. At the baseline operation of 990 rpm, the efficiency is 63.96%, while the Nusselt number is 24.89. If the rotating speed is increased to 1800 rpm, the thermal efficiency decreases to 42.23%. However, the Nusselt number goes to 36.74. In other words, a 50% increase in the convection heat transfer rate can be achieved at a 20% drop in efficiency. When a rotating speed of 1350 rpm is used the efficiency is reduced by 8.6% and allows for an increase in the calculated Nusselt number by 15%. The fan thus plays an important role in enhancing the flow mixing and convective heat transfer in the furnace. Additionally, by controlling the fan speed, it is also possible to improve the temperature uniformity inside the chamber.
Energies 2020, 13, x FOR PEER REVIEW 20 of 23 24.89. If the rotating speed is increased to 1800 rpm, the thermal efficiency decreases to 42.23%. However, the Nusselt number goes to 36.74. In other words, a 50% increase in the convection heat transfer rate can be achieved at a 20% drop in efficiency. When a rotating speed of 1350 rpm is used the efficiency is reduced by 8.6% and allows for an increase in the calculated Nusselt number by 15%. The fan thus plays an important role in enhancing the flow mixing and convective heat transfer in the furnace. Additionally, by controlling the fan speed, it is also possible to improve the temperature uniformity inside the chamber.

Conclusions
A numerical investigation has been carried out, presenting the thermal and fluid dynamics behavior of an electric tempering furnace. A set of conclusions have been obtained on both items and are summarized in what follows.
The simulations predict the thermal homogeneity inside the process chamber effectively.

Conclusions
A numerical investigation has been carried out, presenting the thermal and fluid dynamics behavior of an electric tempering furnace. A set of conclusions have been obtained on both items and are summarized in what follows.
The simulations predict the thermal homogeneity inside the process chamber effectively. Thermal homogeneity was negatively affected by the charging basket use and it was improved by the increase in the rotating speed of the fan. The thermal homogeneity can be associated with the strong recirculation zone formed caused by the rotational force of the fan. The recirculation enhances the flow mixing, increases the residence time of the air in the process chamber, and thus enhances the convective heat transfer.
The calculated air mass flow shows that mass recirculation and the vortex location inside the recirculation zone are affected by the rotating speed variation. When the rotating speed increases more mass is recirculated in the process chamber and the central zone of the recirculation is displaced downward. However, the ratio of mass recirculation to total mass decreases as a result of a higher overall mass flow rate. These values were lower than unity and they were also lower than others found in the literature, suggesting that the fluid dynamics characteristics of the tempering furnace could be further improved.
The thermal efficiency of the tempering furnace could be penalized by an increment in the rotating speed of the fan. This is due to the increase in the fan's mechanical power consumption when the rotating speed increase. Nevertheless, it was shown that thermal efficiency decreases by only 20% when the rpms are doubled, while the heat transfer rate to load is increased by up to 50%. For a furnace operated continuously, this could be translated into lower processing times, and a higher production output, even if that means a lower thermal efficiency operational point. Moreover, it was shown that a higher rotating speed is related to a more homogenous temperature distribution, which is a highly desirable feature in heat treatment equipment.  Acknowledgments: Authors want to thank Forjas Bolivar S.A.S for allowing them to carry out the study.