Three-dimensional fluid-solid coupling heat transfer simulation based on the multireference frame for a side-blown aluminum annealing furnace

In this study, a three-dimensional (3D) numerical simulation model for the flow and heat transfer in a side-blown aluminum annealing furnace (SAAF) is successfully established. Based on the vivid evolutions of the flow field and temperature field, it is confirmed that multiple vortices among the radially distributed nozzles play a key role in reducing the interior flow resistance of the SAAF. The simulated flow distribution agrees remarkably well with the on-site experimental data, which reveals that the model based on the multireference frame method is suitable for describing the 3D fluid-solid coupling heat transfer process inside the SAAF. In addition, to resolve the appreciably uneven temperature distribution in the SAAF, a scheme of guide plate arrangement around the fan is developed to adjust the flow pattern and facilitate the reasonable allocation of the nozzle flow. The standard deviation and the coefficient of variation are obviously declined (∼12%) for both low and standard airspeeds, thereby suggesting that the uniformity of the nozzle flow distribution are expectedly improved. The progress made so far is a substantial step toward achieving high quality, high efficiency and energy savings in aluminum production.


Introduction
With the ongoing increase of the consumption of aluminum products in the construction, aerospace and food packing fields, various production lines of aluminum foils, strips and coils have been produced and implemented. Each stage of aluminum production, such as melting, heating and annealing, consumes a huge amount of energy. As was found in our previous work (Qiu, Feng, Chen, Li, & Zhang, 2018) and others (Carmona & Cortes, 2015;Nieckele, Naccache, & Gomes, 2011;Wang, Zhou, Yan, & Zhou, 2014), the melting process for the regenerative aluminum melting furnace can be accurately described using simulation approaches and a multiobjective optimization method, which is proved to be a good tool for significant energy consumption reductions. Similarly, the annealing process is validated to be useful in improving the mechanical and physical performances of aluminum products (Tsuji, Ito, Saito, & Minamino, 2002). To boost the quality and the efficiency of the aluminum annealing process, a number of researchers have initiated studies on the heat transfer process and established temperature field models. For example, a three-dimensional (3D) model for a continuous CONTACT Yanhui Feng yhfeng@me.ustb.edu.cn annealing furnace is established by Depree et al. (2010) to study the temperature distribution of steel strips. The model based on a short solution time can simultaneously optimize the heat treatment quality, plant throughput and energy consumption. Aiming to elucidate the interlaminar bonding phenomenon of aluminum coils in hydrogen bell-type annealing furnaces, the temperature field and stress distribution of coils are theoretically studied by Saboonchi, Hassanpour, and Abbasi (2008) via establishing a numerical calculation model. Another heat transfer model containing the effects of the heat storage of the workpiece and the heat dissipation from the furnace wall is constructed using the simulation approach by Kang and Rong (2006). Mehta and Sahay (2009) established a heat transfer-phase transformation integrated model of aluminum coil annealing and analyzed the furnace productivity. However, due to the gap between the layers of the aluminum coil, and the oil film formed by the aluminum coil during the rolling process, the contact thermal resistance between the aluminum coil layers is increased, thereby causing the uncertainty of the radial thermal conductivity of the aluminum coil. In recent years, many scholars have established a physical model for heat transfer of aluminum coils, and studied the equivalent radial thermal conductivity of aluminum coils. A reasonable and effective heat transfer model based on a large amount of on-site production data was proposed by Chen and Gu (2007) to predict the temperature field of a typical annealing furnace for aluminum coils production. Mathematical expressions for the equivalent radial thermal conductivity (κ e,a ) of the aluminum coils and the total gas circulation flow are obtained based on regression analysis. The reason why hydrogen is superior to nitrogen in the annealing process is profoundly studied by Pal, Datta, and Sahay (2005), and a theoretical formula for κ e,a is deduced by using the network diagram method.
However, most of the above studies are based on traditional radial heating furnaces. The heat transfer efficiency of the traditional radial heating annealing furnace is estimated to be 60% according to a layer-by-layer heat transfer mechanism initiated by the heat convection from the hot air to the outermost layer of the coil. Due to the gap and oil film formed between the coil layers in the rolling process, the κ e,a value of the coil is extremely small, which seriously hinders the heat transfer efficiency of the annealing process. In comparison to traditional radial heating annealing furnaces, the temperature distribution within the aluminum coils is more uniform in the same boundary conditions by adopting the side-blown heating method. Hence, a novel side-blown aluminum annealing furnace (SAAF) is proposed in this paper. In addition, the 3D fluid-solid coupling model based on the multireference frame (MRF) is established to simulate the annealing process of aluminum coils, aiming to improve the annealing quality and efficiency.

Physical model
As shown in Figure 1, the SAAF is mainly composed of a circulating fan located on the top, radiation tubes and nozzles on both sides of the furnace, and horizontal and vertical baffles around it. The reel refers to the central axis of the aluminum coil, and the furnace is a 3D space surrounded by baffles in the furnace for annealing and heating of aluminum coils. The furnace structure of the annealing furnace is described using the Cartesian coordinate system. The origin of the coordinates is at the center of the aluminum coil, and the x-axis coincides with the aluminum reel line and points to the right nozzle. The y-axis is perpendicular to the aluminum reel line and points to the fan located in the upper part of the furnace. In this paper, the Integrated Computer Engineering and Manufacturing code for Computational Fluid Dynamics (ICEM-CFD) is employed for modeling and meshing. In this paper, the whole model is solved, and the whole calculation area is divided into two parts, the rotating zone and the stationary zone. Due to the complexity of the fluid region model, the Tet/Hybrid element and unstructured mesh are used to divide the grids. Particularly, the prism element near the wall is applied to satisfy the boundary layer conditions. What's more, all grids are coupled at the interface to form integrated model meshes, and the total number of grids is 5,031,584. The grid model is imported into ANSYS-CFX and then loaded for calculation.

Mathematical model
In recent years, numerical computing technology has been continuously developed and applied to solve many complicated problems (Akbarian et al., 2018;Ghalandari, Koohshahi, Mohamadian, Shamshirband, & Chau, 2019;Mosavi, Shamshirband, Salwana, Chau, & Tah, 2019;Mou, He, Zhao, & Chau, 2017). Here, a 3D physical model for the SAAF is established, and an MRF model is used to increase the accuracy of the calculation. When the MRF is used, the entire computing domain is divided into three subdomains, that is, the rotating impeller zone (rotating zone), aluminum coil region and the rest of the stationary zone. Each subdomain has its own motion forms: stationary, rotating or translational. Whether it is side spray heating or radial heating, the division of the calculation domain will not be affected. The governing equations of the flow field are solved in each subdomain, and the flow field information is exchanged by converting the relative velocity into the absolute velocity at the interface of the subdomain. The realizable k-turbulence model and the wall function are used to describe the flow field.
To simplify the calculation, the following hypothesis were made on the numerical model: (1) neglecting the radiative heat transfer due to the small emissivity of aluminum and high flow rate of air, leading to the fact that the heat transfer in the SAAF relies mainly on convective heat transfer; (2) the thicknesses of both the baffles and the fan blades are small enough to be treated as a nothickness wall; and (3) assuming that the outer wall is adiabatic, which means that no heat is transferred in or out of the system.
The governing equations are 3D Navier-Stokes equations, including mass conservation, momentum conservation and energy conservation equations. They are discretized by using the finite volume method based on the finite element method (FEM), and thus, the method retains the numerical accuracy of FEM and the conservation properties of the finite volume method (Ramezanizadeh, Alhuyi Nazari, Ahmadi, & Chau, 2018).

Multireference frame model
The MRF model is the simplest of the three approaches for multiple zones, i.e. the mixing plane model, the multireference frame model and the sliding mesh model. It makes steady-state approximation calculation for each element with different rotating or moving speeds. In addition, the MRF model is used to deal with the mixing and rotation problems, which can effectively improve the computational accuracy (Bujalski, Jaworski, & Nienow, 2002;Chapman, Sudhoff, & Whitcomb, 1998;Cruz & Cardoso, 2005). When running the MRF model, the zero equation turbulence model will not suffice since it does not properly capture the swirl effects of the rotor blades. It is necessary to adopt the realizable two-equation turbulence model.

Realizable k-turbulence model
The realizable k-model is adopted in the turbulence model, which is widely used in engineering and other applied sciences (Hargreaves & Wright, 2007;Seeta Ratnam & Vengadesan, 2008). Compared with the standard k-turbulence model, the main difference of the realizable k-model is the increase in the calculation of turbulent viscosity and the transfer equation that changes the diffusion rate. In the realizable k-model, the dynamic viscosity calculation coefficient is no longer constant, but is related to the strain rate. In addition, the calculation of the turbulent dissipation rate is also different from the standard k-model. At present, the model has been widely used to calculate more complicated flow conditions such as a swirling shear flow, a planar mixed flow (Bridgeman, Jefferson, & Parsons, 2009), a jet (Li, Wang, & Chen, 2018;Zhang, Yang, & Li, 2019), and so on: where c μ represents the dynamic viscosity coefficient. The coefficients U* and c 1 are both related to the timeaveraged strain rate, where U* is also related to the timeaveraged rotation rate and used to describe the effect of the rotating flow field.

Wall function method
In the near-wall region, the flow is not fully developed by the wall limitation, so the realizable k-model is no longer suitable for the calculation because of the insufficient development of the turbulence. To solve the flow calculation in the near-wall region, a wall function method is adopted. It is believed that the abovementioned method is the most popular way to account for wall effects (Knopp, Alrutz, & Schwamborn, 2006). It can be expressed by the following expressions: where u + , y + and T + are three nondimensional parameters that denote the speed, distance and temperature, respectively. K p and T p are the turbulent kinetic energy and temperature of the node P, respectively. y p represents the distance from node P to the wall and μ represents the kinetic viscosity of the fluid. E is the experimental constant, and E = 9.8 in the present simulation. T w and q w are the temperature and heat flux on the wall, respectively.

Radiation tubes
The radiation tube is also called immersion heating tube. The working principle is that air sweeps through the high-temperature radiation tube wall and then heated. In this process, since air does not have the ability to emit and absorb radiant energy within the annealing temperature range of the aluminum coil, the radiative heat exchange between the air and the radiant tube is not considered, and the heat transfer is mainly in the form of convective heat transfer. It is widely used in heating furnace and is the main component of heating furnace. In this paper, the condition of constant wall temperatures is adopted for the radiation tubes, which is divided into two parts. The temperature of the upper part is set to 1200 K, while the lower part is 1000 K.

Fan
The air enters the furnace through the centrifugal fan. And the whole calculation area includes the rotating zone and the stationary zone. The MRF model is used to couple the rotating vane with the static wall. Figure 2 shows the meshing results of the fan. The velocity inlet boundary condition and the pressure outlet boundary condition are applied to the fan inlet and fan outlet, respectively. In this paper, the rated speed of the fan used in the SAAF simulation is 216 rad/s and the rated flow is 1,370,000 m 3 /h. And the aluminum coil is annealed under an air atmosphere.

Wall
The heat transfer mode between the inner surface of the SAAF and the air is dominated by heat convection, and the coefficient of convective heat transfer varies with the temperature of the air in the furnace. The symmetrical boundary conditions are applied to the walls, which are parallel to the axial of the fan. In addition, the adiabatic boundary condition is applied to the outer wall.

Aluminum coil
A mathematical model is developed based on the theoretical analysis for the description of the aluminum coil and the corresponding computational model is built via the ICEM-CFD software. ANSYS 16.0 is used to numerically solve the mathematical model of annealing furnace, and a separate solver is used to calculate the steady state. Here, the heat transfer of the aluminum coil is calculated, which does not involve any flow calculation. The physical parameters of the aluminum coil are shown in Table 1.
Due to the gap between the aluminum coil and the oil film, the uncertainty of the radial thermal conductivity of the aluminum coil arises spontaneously. In this paper, the radial thermal conductivity of the aluminum coil is calculated by the formula proposed by Terrola (1989). The multi-reference system model is used to divide the entire calculation area into multiple different areas. The polar coordinate reference system is used in the area where the aluminum coil nozzle is located. The axial thermal conductivity of the aluminum coil is a fixed value, and the C++ self-programming is used to define the radial thermal conductivity of the aluminum coil which varies with temperature. The continuity assumption is made for the aluminum coil, and the two-dimensional diagram of the aluminum coil is shown in Figure 3. The heat transfer of the aluminum coil can be calculated according to the following formula: where c p is the specific heat capacity of the aluminum coil. λ x and λ r represent the axial and the radial thermal conductivity of the aluminum coil, respectively.

Evolution of the flow field
The velocity vector of the z = 0 cross-section is shown in Figure 4. It is clearly seen that the gas enters the fan  nozzles. Due to the low velocity of the jet from the nozzles, there is no strong impact on the end of the aluminum coil. All nozzles are distributed on the ring surface parallel to the aluminum coil ends, and no nozzles are facing the reel, which makes the air not flow along the axial direction of the aluminum coil. It is found that the flow rate of the radial rotating flow in the reel is very low, which can be reasonably neglected.
It is clearly seen that there are multiple vortices in the flow field from the local magnification images. Figure 5(a) shows the local magnification representation of the velocity vector at position 1 as labeled in Figure 4. Because the outlet velocity of the fan is relatively high, there is air entrainment and eddy formation at the lower part of the fan. The flow field in the inner ring nozzles is demonstrated in Figure 5(b). Because of the change of the flow direction near the nozzle, the particles of the fluid collide with each other and thus produce vortices. The collision process is accompanied by severe friction and momentum exchange processes, resulting in energy loss, which further hinders fluid flow, such as the reverse vortex between inner ring nozzles. Of note is the fact that the vortices hamper the flow and thus do not promote the heat transfer in the SAAF. Moreover, multiple vortices occur between the nozzles located at different radial locations due to the jet effect ( Figure 5(c,d)), and the existence  of vortices increases the disturbance of the fluid, enabling the improvement in the convective heat transfer between the airflow and the aluminum coil.

Evolution of the temperature field
The temperature distribution of the z = 0 and x = 0 cross-sections are shown in Figure 6(a). According to the picture, we can conclude that the temperature distribution in the SAAF is not uniform. Specifically, the temperature at the bottom of the SAAF is relatively low, while the temperature inside the coil is much lower, due to the internal gas being almost stationary and the poor heat exchange between air and the aluminum coil. The difference between the jet temperature of nozzles and the low-temperature zone of the reel center is approximately 100 K. The heating of the aluminum coil mainly relies on the convective heat transfer of the end faces and the outer circumference.
The high-and low-temperature areas of the aluminum coil are clearly shown in Figure 6(b). The hightemperature area is distributed on the outer peripheral surface away from the circumference of the reel, while the low-temperature area appears on the inner circumference near the reel, which is the direct result of no gas flow in the internal part of the reel. Moreover, because of the suction of the fan, the top of the SAAF renders a higher wind speed than the bottom, and thus the upper part of the aluminum coil is heated much faster.
To monitor the temperature variation of the aluminum coil with time, four positions are selected as temperature monitoring points for the aluminum coil. As shown in Figure 7, P 1 , P 2 and P 3 are the monitoring points of the right end of the aluminum coil and P 4 is the monitoring point of the middle plane. The temperature evolutions of the four monitoring points with time are summarized in Figure 8. For monitoring points at different locations, the temperature will definitely vary with the annealing time. The temperature change of the monitoring point on the end face of the aluminum material is consistent. Hence, the temperature changes of P 1 , P 2 , and P 3 are basically the same. However, since the monitoring point P 2 is located at the inner nozzle, the heating rate is slightly slower than that of P 1 and P 3 . The annealing furnace studied in this paper uses side spray heating to heat the aluminum coil, and the temperature on the end Figure 7. Schematic diagram of monitoring points for the aluminum coil (P 1 , P 2 and P 3 are the monitoring points on the right end of the aluminum coil and P 4 is the monitoring point on the middle surface). face of the aluminum coil is first raised and then transmitted to the inside. Therefore, the temperature at the center section of the aluminum coil is lower than the end surface temperature. Since the radiation tubes are treated with the first kind of boundary condition, the aluminum coil is heated continuously during 12 h for the numerical simulation. During the first half hour, the gas temperature near the outlet of the fan rises rapidly, which is followed by a slow climb. With the increasing heating time, the temperature differences between the monitoring points of the right end (hot point) and the center plane (cold point) decrease gradually.

Optimization and analysis of structural parameters: the guide plate
Note that the above-mentioned physical model cannot be verified with the actual temperature field data since the guide plate is not included in the present model while it is used in actual apparatuses. It is necessary to design a flow diversion structure in the upper space of the SAAF to weaken the influence of the vortex generated by the fan blade rotation. The optimization result in this section would be an effective guide for reducing the power dissipation of the fan.

Simulation and verification of the flow field with the empty furnace
As verification, the velocity distribution in an empty furnace with a guide plate installed around the fan has been calculated. Figure 9(a,b) illustrate the velocity distribution on the z = 0 and x = 0 cross-sections in the empty SAAF, respectively. The cross-section of the fan with the guide plate installation in the SAAF is shown in Figure 9(c). Clearly, the entrance speed is the highest due to the fan. The jet from the nozzles also causes a high velocity nearby, which makes the verified results rational and reliable.
After a certain number of iterations, the air quantity due to the fan suction generally stabilizes (Figure 10(a)). Although a slight fluctuation in the velocity occurs in the central axes of nozzles 1 and 2, the flow field can be considered to be nearly stable (Figure 10(b,c)).
The article compares the center velocity data of the field test with the simulation results at room temperature to verify the correctness of the simulation method and the reliability of the results. In the case of an empty furnace and a stable operation of the fan, a hot spot anemometer is held and placed close to the spout to measure the velocity of the center of the spout. The measured values and the numerical simulation results of the velocity along the central axes of the nozzles are summarized in Table 2. To make the computational results more visible, a comparison of the graphs of the center axis velocities of the inner, middle and outer nozzles is shown in Figure 11. It can be concluded from Table 2 and Figure 11 that the center axis velocities of the inner ring, middle ring and outer ring nozzles are close to each other.
The relative deviation between the simulation values and the experimental values of each point (E d ) is calculated based on the statistical method, which is given by where v s and v t denote the simulation results and the experimental values of the center axis velocity for each nozzle, respectively. The average relative deviation (E d ) of the center axis velocity of each nozzle and the standard deviation (S d ) are calculated using the following two expressions, respectively: The comparison of the simulation velocity and experimental values from the corresponding positions of the empty furnace indicates that the relative deviation is small, which well validates the reliability of the simulation. Thus, it could be used to analyze the issues regarding the heat transfer and airflow.

Flow field in the furnace under loading at low wind speed
In this section, the flow field in the annealing furnace and the flow distribution of various nozzles at different  radial positions are investigated at a low speed of approximately 66.4 rad/s. The velocity nephograms of the fan cross section, the nozzle section and the cross section of the furnace in different heights with and without the guide plate are obtained. As shown in Figure 12, the air velocity in the area between the furnace top and the circulating fan is higher. Another high-speed zone is located between the two sides of the furnace body and the guide plate. Furthermore, the existence of the guide plate increases the air velocity near the radiation tube area, which leads to the enhancement of the convective heat transfer effect between the radiation tube and the air.
The velocity cloud diagram of the cross sections of the furnace in different heights under two structures, that is, with and without the guide plate, is summarized   in Figure 13. It can be clearly seen that the maximum velocity occurs at the nozzle of all height directions for both structures. The velocity cloud diagrams at the height of y = 0 m clearly show the flow field inside the nozzles at different radial positions. The cross sections at the height of y = 0.38 m and y = −0.38 m can only reflect the flow field inside part of the nozzles. Figure 14 demonstrates the velocity distributions of various nozzles with or without guide plates. It can be clearly seen that the velocity distribution with guide plates is more uniform than the case without the guide plate installation. By detecting the variation of fan inlet air volume with iteration steps, it can be obtained that the fan inlet volume flow before and after the installation of the guide plate is 11.6 and 11.3 m/s, respectively. Therefore, the existence of the guide plate will make the fan inlet velocity decrease slightly. As a result, the maximum velocity of the nozzle decreased after installation of the guide plate. The underlying mechanism for this phenomenon can be further explained by the flow distribution ( Figure 15). And it can be seen from Figure 15 that before the addition of the guide plate, the flow distribution of the nozzles at different radial positions differs greatly, which is more pronounced for the nozzles at the upper left.
Clearly, the flow rate of the bottom nozzles is less than the upper counterpart because of the guide plate, which validates that the guide plate plays a key role in developing the uniform flow distribution of nozzles along the radial direction. It is necessary to compare the discreteness of two sets of data. If the units or the average values are different, S d cannot be used, and instead the ratio of the standard deviation to the average, i.e. the coefficient of variation (C.V), is used. The equation for the C.V is defined by The S d and C.V values of the nozzle flow at different radial locations are compared with or without guide plates, as shown in Table 3.
It is clearly seen that the C.V values of both the middle and outer ring nozzles appreciably decrease after the guide plate is used. In contrast, the C.V of the inner ring nozzle increases by 23.54%. It is concluded that theexistence of the guide plate makes the flow distribution of the nozzles in the middle and outer ring more uniform, while the uniformity in the inner ring decreases slightly. Note that the overall nozzle flow distribution becomes more uniform.

Flow field in furnace under loading at standardized wind speed
The flow distributions of the nozzles with and without guide plates at the standardized airspeed ( Figure 16) are coincident with that at a low speed. Generally, the existence of the guide plate improves the uniformity of the flow distribution of the nozzle. The S d and C.V values with or without guide plates at the standardized airspeed are summarized in Table 4. Similar conclusions are obtained, which further verifies the correctness of the results. The standard deviation and variation of coefficient with the guide plate are compared under two conditions (low wind speed and standardized wind speed) as shown in Figure 17. It can be clearly seen that both the standard deviation and the variation of coefficient under the two conditions declined after installing the guide plate. The standard deviation and the coefficient of variation reflect the uniformity of the flow distribution. Therefore, it is safe to conclude that the flow distribution of all the nozzles is more uniform than that of the nozzles without the guide plate.

Conclusions
In this paper, based on the multi-reference system, a numerical model for the coupling of air-aluminum coil flow and heat transfer is established. The numerical simulation and structural optimization analysis of the annealing process of the aluminum coil in the SAAF are carried out. The airflow and heat transfer in the side-blown aluminum annealing furnace are analyzed in Figure 16. Flow distribution of various nozzles at the standardized airspeed.
detail, and the evolution characteristics of the flow field and temperature field are obtained. Moreover, the role of the guide plate on the nozzle flow distribution is investigated. By comparing the simulation results with the field data, the relative deviation and coefficient of variation are used to characterize the uniformity. The main conclusions are listed as follows.
(1) For the side-blown aluminum annealing furnace, the reverse vortices formed among the inner ring nozzles increases the flow resistance. However, the other vortices formed among the nozzles at different radial positions enhance the heat transfer between the fluid and the aluminum coil.
(2) The guide plate arranged around the fan at a low airspeed enables the flow distribution of the nozzles in the middle and outer rings to be more uniform, while the uniformity of the inner ring region decreases slightly, which is clearly reflected by the standard deviation and the coefficient of variation. (3) The guide plate also improves the uniformity of the flow distribution of the nozzles under the standard airspeed, which further demonstrates the validity of the simulation results. (4) There is almost no flow in the lower part of the furnace, and thus the heating speed of the lower part of the aluminum coil is obviously slower than that of the upper part. Therefore, it is recommended to combine the traditional radial heating with the side spray heating. Furthermore, to achieve fast and even heating, the intake air volume of the two heating modes needs to be optimized according to the optimization design. (5) This manuscript only studies the annealing process of aluminum coils in side spray annealing furnaces without guide plates. It is desirable to further study the uniformity of the temperature of the aluminum coil before and after the installation of the baffle in the subsequent researches. Moreover, the factors affecting the annealing process should be discussed and analyzed, and the research on system power consumption will be carried out in future work.