**1. 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, amongs<sup>t</sup> 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 usually achieved by heating the steel right after the heat treatment, following a controlled temperature

rise to a preset level in the 670–970 K range [6]. If the temperature is not uniform during the heating processes, some parts of the heat-treated steel might have an undesirable structure and mechanical properties. This might be unacceptable for highly demanding applications—in the aeronautics industry for example.

Temperature uniformity inside heat treatment furnaces is commonly rated according to the DIN 17052-1 or the AMS2750 standards [7,8]. Both standards define "furnace classes" based on the maximum temperature deviation within the workspace of the furnace concerning the temperature setpoint. The most stringent class only allows for a ±3 K temperature di fference in such a uniformity range. This range is determined by following a time-consuming and expensive experimental survey that only measures temperature by a few points inside the furnace.

Convective furnaces use fans to recirculate a hot fluid and to achieve a homogeneous temperature. Previous research has shown that such furnaces require a preheating stage where energy consumption is the highest and it is proportional to the preheating time [1]. The results of that study also showed that using ba ffles or other obstacles inside the furnace help reduce the preheating time.

Many studies concerning the optimization of furnaces have been oriented to maximize the quality of the product while minimizing the amount of energy required. In the food industry, for example, the quality of cooking can be controlled by keeping the uniformity of the temperature in the oven [9]. In this case, it was shown that the heating coil placement on a side position guaranteed the homogeneity of the temperature inside the furnace. Reference [10] reported that the uniformity of temperature in a laboratory drying oven was improved experimentally through changes in the ventilation system. By using CFD simulations, the influence of the position of the heaters and deflectors to redirect the air-flow more e ffectively was analyzed. The heated air was distributed directly on the device's storage space and the authors analyzed the e ffect of fan speed on the temperature distribution. Four di fferent values were tested, i.e., 1500 rpm, 2000 rpm (existing fan), 2500 rpm, and 3000 rpm. The results led to the conclusion that the higher the rotating speed of the fan, the lower the obtained temperature di fference. However, the achieved improvement was only about 0.1 K, which can be considered to be a very low value. These studies recommended the use of CFD to create geometric prototypes that guarantee the uniformity of temperature in convection ovens.

Simulations of a furnace chamber during heat treatment using natural gas combustion have been reported in the literature, where a boundary condition model for the estimation of heat flux through the walls was proposed and validated [11]. However, a finite-element method, rather than a finite-volume method, was used, and only the energy transport equation was solved. In another study [12], an electric nitriding furnace was simulated using a commercial CFD package and velocity uniformity was evaluated using a spatial criterion based on the mean and actual values of the velocity magnitude. There is also another simulation study regarding a nitriding furnace, where a conjugated heat transfer approach was used [13].

In this research, the temperature distribution, fluid dynamics, and heat transfer behavior in an electric tempering furnace are studied. Thus far we know that there is limited research concerning these characteristics in this specific kind of furnace. Numerical simulations have conducted using a relevant design factor as the rotating speeds of the fan. In addition, thermal e fficiency has been calculated when the rotating speed is changed. Taking into account the fact that internal velocities and temperatures fields inside the furnace are di fficult to directly measure, the numerical approach allows for the identification of the main phenomena governing the thermal and fluid behavior and could contribute to the furnace design for process improvements.

#### **2. Numerical Simulations**

#### *2.1. System Description*

The study investigates a top loading furnace commonly referred to as a pit-type furnace (Figure 1). Pit furnaces have been widely used for steel tempering (394 K to 1033 K) due to the easy loading and 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.

**Figure 1.** Tempering furnace studied in this work. (**a**) Isometric external view of the furnace and (**b**) cut view of the furnace with the coordinate system.

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 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.

**Figure 2.** Setpoint temperature at the cooling (120–150 min) and sustain stage (150–390 min) of the tempering process.

#### *2.2. 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 lowMach 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

$$\frac{\partial \overline{\rho}}{\partial t} + \nabla \cdot (\overline{\rho} \overline{v}) = 0,\tag{1}$$

where ρ and *v* are the average density and velocity. The momentum equation is also solved:

$$\frac{\partial}{\partial t}(\widetilde{\rho v}) + \nabla \cdot (\widetilde{\rho v} \widetilde{v}) = \nabla \cdot \left\{ \mu\_l \left[ \left( \nabla \widetilde{v} + \nabla \widetilde{v}^T \right) - \frac{2}{3} \nabla \cdot \widetilde{v} l \right] \right\} - \nabla \cdot (\widetilde{\rho v} \widetilde{v} \widetilde{v}) - \nabla \widetilde{p} + \widetilde{\rho} \widetilde{g}, \tag{2}$$

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*<sup>2</sup> *v* is a term that must be modeled. In this work, the Boussinesq hypothesis was followed:

$$-\overline{\rho}\widetilde{\upsilon\upsilon}^{\prime}\overline{\upsilon}^{\prime} = \mu\_{l}(\nabla\overline{\upsilon} + \nabla\overline{\upsilon}^{\prime}) - \frac{2}{3}(\overline{\rho}\overline{k} + \mu\_{l}\nabla\cdot\overline{\upsilon}l),\tag{3}$$

where *k* is the averaged turbulent kinetic energy:

$$\widetilde{k} = \frac{1}{2} \Big[ \widehat{\left(u\right)}^2 + \widehat{\left(v\right)}^2 + \widehat{\left(w\right)}^2 \Big] = \frac{1}{2} \sum\_{i=1}^3 \widehat{\left(v'\_i\right)}^2. \tag{4}$$

The turbulent viscosity μ*t* is obtained from the following equation:

$$
\mu\_t = \overline{\rho} \mathbf{C}\_{\mu} \frac{\overline{k^2}}{\overline{\varepsilon}},
\tag{5}
$$

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 *Energies* **2020**, *13*, 3655

successfully in the past for simulations of thermal treatment furnaces [16]. In this model, two additional transport equations for *k* and ε must be solved:

$$\frac{\partial}{\partial t} \left( \overline{\rho} \overline{k} \right) + \nabla \cdot \left( \overline{\rho} \overline{v} \overline{k} \right) = \nabla \cdot \left[ \left( \mu\_l + \frac{\mu\_l}{\sigma\_k} \right) \nabla \overline{k} \right] - \overline{\rho} \overline{v'v'} : \nabla \overline{v} - \overline{\rho} \overline{\varepsilon}\_{\prime} \tag{6}$$

$$\frac{\partial}{\partial t}(\overline{\rho}\overline{\varepsilon}) + \nabla \cdot (\overline{\rho}\overline{\upsilon}\,\overline{\varepsilon}) = \nabla \cdot \left[ (\mu\_l + \frac{\mu\_l}{\sigma\_\varepsilon}) \nabla \overline{\varepsilon} \right] - \mathbb{C}\_{\varepsilon1} \overline{\overline{\rho}} \, \frac{\overline{\varepsilon}}{\overline{k}} \overline{\overline{\upsilon}\prime\prime} : \nabla \overline{\upsilon} - \mathbb{C}\_{\varepsilon2} \overline{\overline{\rho}} \, \frac{\overline{\varepsilon}^2}{\overline{k}} \, \tag{7}$$

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:

$$
\frac{
\partial
}{
\partial t
}
\left(
\overline{\rho}\overline{E}
\right) + \nabla \cdot \left(
\overline{v}
\left(
\overline{\rho}\overline{E} + \overline{p}
\right)
\right) = \nabla \cdot \left(
k\_l \nabla \overline{T}
\right) + S\_{rad\,\prime}
\tag{8}
$$

where *kt* is the thermal conductivity of the air, *T* is the average temperature and *Srad* is the source term due to radiation. The total non-chemical energy E is defined as

$$E = h - \frac{p}{\rho} + \frac{1}{2}v^2,\tag{9}$$

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*):

$$\frac{\partial}{\partial s}l\left(\stackrel{\rightharpoonup}{s}\right) + (a + \sigma\_5)l\left(\stackrel{\rightharpoonup}{s}\right) = an^2\sigma \frac{T^4}{\pi} + \frac{\sigma\_s}{4\pi} \int\_0^{4\pi} l\left(\stackrel{\rightharpoonup}{s}\right) \phi\left(\stackrel{\rightharpoonup}{s}, \stackrel{\rightharpoonup}{s}'\right) d\Omega,\tag{10}$$

where *<sup>I</sup>*<sup>→</sup>*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−<sup>8</sup> Wm−<sup>2</sup> <sup>K</sup>−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]:

$$S\_{rad} = aG - 4an^2 \sigma T^4,\tag{11}$$

where *n* is the refractive index and *G* is the incident radiation, which is evaluated as

$$G = \int\_{4\pi} Id \Omega. \tag{12}$$

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.

**Figure 3.** Fan modeled in this work. (**a**) Isometric view; (**b**) cylindrical boundary and domain configured with rotating speed using a moving-reference frame (MRF).

#### *2.3. 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.


**Table 1.** Wall boundary conditions for baseline operation.

In this work, the heat flux .*q* was assumed to be due to a combination of convection and radiation:

$$
\dot{q} = \dot{q}\_{conv} + \dot{q}\_{rad}.\tag{13}
$$

The tempering furnace studied in this work is located inside a room with no significant airflow coming in or out. Therefore, a 5 Wm−<sup>2</sup> K−<sup>1</sup> natural convective coefficient *hconv* was assumed [23], and the convective heat flux .*qconv* was calculated with Newton's law of cooling:

$$
\dot{q}\_{conv} = h\_{conv}(T\_s - T\_{surr}),
\tag{14}
$$

where the surrounding temperature *Tsurr* is 300 K on average and *Ts* 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]:

$$
\dot{q}\_{\rm rnd} = \varepsilon \sigma (T\_s^4 - T\_{surr}^4)\_\prime \tag{15}
$$

where ε is the surface emissivity. All surfaces were considered opaque with a di ffusive 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−<sup>6</sup> 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.

**Figure 4.** Points selected for temperature measurement on the exterior of the furnace walls (measurements in m).

#### *2.4. 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.

**Figure 5.** Grid systems (**a**) coarse mesh, 1.5 million cells, and (**b**) fine mesh, 3 million cells.

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.

#### *2.5. Convective Coe*ffi*cient 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/m3, and specific heat of 477 J/(kg K) [23]. Calculated from density and volume, the mass of this basket is 54.7 kg.

**Figure 7.** Stainless steel charging basket. (**a**) Empty basket; (**b**) loaded with parts to be tempered.

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/m3, 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]:

$$N\mu = \mathbb{C}\,R\varepsilon^m Pr^{1/3}\,,\tag{16}$$

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]:

$$Nu = \, \mathbb{C} \, Re^m Pr^n \Big(\frac{Pr}{Pr\_s}\Big)^{1/4},\tag{17}$$

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:

$$Nu = 0.3 + \frac{0.62 Re^{1/2} \text{Pr}^{1/3}}{\left[1 + \left(0.4/\text{Pr}\right)^{2/3}\right]^{1/4}} \left[1 + \left(\frac{Re}{282,000}\right)^{5/8}\right]^{4/5}.\tag{18}$$

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

$$\eta = \frac{\text{Heat transfer rate to the load}}{\text{Electric power input}},\tag{19}$$

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]:

$$
\left(\frac{\dot{\mathcal{Q}}}{\omega D^3}\right)\_1 = \left(\frac{\dot{\mathcal{Q}}}{\omega D^3}\right)\_2 = \Phi,\tag{20}
$$

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 *CP*, then the following similarity law also should apply:

$$
\left(\frac{\dot{W}\_{shaft}}{\rho \omega^3 D^5}\right)\_1 = \left(\frac{\dot{W}\_{shaft}}{\rho \omega^3 D^5}\right)\_2 = \text{C}\_{P\prime} \tag{21}
$$

where *Wsha 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.

#### **3. 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.

#### *3.1. 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 agreemen<sup>t</sup> 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.

**Figure 8.** Temperature distribution of the baseline operation without the charging basket. (**a**) Cut view of the symmetry plane; (**b**) cut views at different heights.

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 augmen<sup>t</sup> 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.

**Figure 9.** Streamlines inside the furnace for baseline operation and air mass recirculation behavior. Velocity legend applies to the recirculating gases.

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 10.** Distribution of Y+ on the walls for the baseline operation without the charging basket. The Y+ legend applies to the gases.

**Figure 11.** Temperature distribution of the baseline operation with a charging basket.

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.

**Figure 12.** Streamlines inside the furnace for baseline operation with a charging basket. The velocity streamline legend applies to the recirculating gases.

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.

**Figure 13.** Sankey diagram for energy flows during the heating stage of the tempering process for the baseline operation of 990 rpm. (**a**) During the preheating stage; (**b**) at steady state.

If the rotating speed is increased from 990 rpm, fan energy demand would also increase, and e fficiency 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 a ffected 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.

#### *3.2. E*ff*ect 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 di fferent 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.

**Figure 14.** Temperature profile at a horizontal line located in the middle of the workspace for the five di fferent angular rotation rates studied.

From Figure 14, the maximum temperature di fference of the baseline operation of 990 rpm, between radial distances of −0.25 and 0.25 m, is 47.5 K. This di fference 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 di fferent 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 di fferentials become larger. As shown before, for the 1800 rpm case, temperature homogeneity is enhanced, and the velocity profiles sugges<sup>t</sup> 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 e ffect in the convective heat transfer.

**Figure 15.** Axial velocity profile at a horizontal line located in the middle of the workspace for the five different rotating speeds studied.

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.

**Figure 16.** Mass recirculation (**a**) and temperature CU (**b**) at several heights 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 di fference at 990 rpm case was 47.5 K and the di fference 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 su ffers 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 (mrmax) 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–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.

**Figure 17.** Relative mass recirculation at a di fferent rotating speed of the fan.

#### *3.3. Thermal E*ffi*ciency and Heat Transfer to the Load*

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 coe fficient Φ for any two di fferent operational conditions. This validates that the fan used in this study follows the classic laws for fans with low uncertainties.

**Figure 18.** Recirculation flow rate of hot air at di fferent rotating speeds of the fan.

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 coe fficient *CP* for any two di fferent 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:

$$
\left(\frac{\mathcal{W}\_{shaft}}{\alpha^3}\right)\_1 = \left(\frac{\mathcal{W}\_{shaft}}{\alpha^3}\right)\_2\tag{22}
$$

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 a ffect the thermal e fficiency of the furnace; however, it should be contrasted with its positives e ffect 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 e fficiency η as a function of the fan rotating speed. Both can be approximated to a quadratic equation as a function of the rpms with a coe fficient of determination R<sup>2</sup> of 0.99. At the baseline operation of 990 rpm, the e fficiency is 63.96%, while the Nusselt number is 24.89. If the rotating speed is increased to 1800 rpm, the thermal e fficiency 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 e fficiency. When a rotating speed of 1350 rpm is used the e fficiency 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.

**Figure 19.** Nusselt number and e fficiency variation with the rotating speed of the fan. N = rpm.

## **4. 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 e ffectively. Thermal homogeneity was negatively a ffected 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 a ffected 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 e fficiency 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 e fficiency 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 e fficiency 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.

**Author Contributions:** Investigation, I.D.P.-C., P.N.A.-T. and L.F.C.-S.; writing—original draft, I.D.P.-C., P.N.A.-T., and L.F.C.-S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** Authors want to thank Forjas Bolivar S.A.S for allowing them to carry out the study.

**Conflicts of Interest:** The authors declare no conflict of interest.
