- Email: [email protected]
Prediction of the heating characteristic of billets in a walking hearth type reheating furnace using CFD
International Journal of Heat and Mass Transfer 92 (2016) 675–688
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Prediction of the heating characteristic of billets in a walking hearth type reheating furnace using CFD Rene Prieler a,⇑, Bernhard Mayr a, Martin Demuth b, Burkhardt Holleis b, Christoph Hochenauer a a b
Institute of Thermal Engineering, Graz University of Technology, Inffeldgasse 25/B, A-8010 Graz, Austria Messer Austria GmbH – Kompetenzzentrum Metallurgie, Industriestraße 5, A-2352 Gumpoldskirchen, Austria
a r t i c l e
i n f o
Article history: Received 17 June 2015 Received in revised form 3 August 2015 Accepted 4 August 2015 Available online 26 September 2015 Keywords: Combustion modelling Reheating furnace Computational fluid dynamics Heating characteristic Steady flamelet
a b s t r a c t A natural-gas fired walking hearth type furnace for the reheating of steel billets was investigated in this work. A novel numerical approach was used to predict the gas phase combustion, heat transfer and transient heating characteristics of the billets in the furnace. In contrast to conventional coupled simulation of the heat transfer and the heating characteristics, the iterative solution procedure used considers steady-state gas phase combustion of the furnace and unsteady simulation of the billets separately. Due to the steady flamelet model for combustion modelling the calculation time was kept to a minimum although a detailed reaction mechanism with 17 species and 25 reversible reactions was used. The reaction mechanism is applicable for different combustion environments and can be used for future investigations of the furnace after adaption for oxygen enriched combustion. Burner geometry was modelled in detail for a high accuracy prediction of the flame shape of the flat flame burners as well as temperature and species concentrations in the vicinity of the burner. The simulation revealed that 93% of the total heat flux to the billets was contributed by radiation. Transient heating showed a maximum temperature difference inside the billets in the heating zone at 192 K decreasing to a value between 42 and 49 K at the end of the furnace. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Large scale reheating furnaces are employed in the steel industry for heating up steel billets or slabs to a desired temperature level for further treatment in a rolling mill. Furthermore, a uniform temperature distribution is desired inside the billets with optimised operating conditions regarding fuel consumption, residence time, NOx and CO2 emissions etc. According to an IPCC report, about 60% of global CO2 emissions are produced by fossil fuels fired processes [1]. In the world today environmental issues and cost efficient production are the subject of major attention in numerous industrial sectors. Significant efforts are being made to increase efficiency against this background. In the past few years as a result of increasing computer power, many researchers have resorted to numerical methods for investigating reheating furnaces as these are applied in the steel manufacturing process. Numerical models can offer a detailed view of the transport phenomena in the furnaces and can be a helpful tool to optimise the heating process, transportation of the furnace loads and other phenomena. ⇑ Corresponding author. Tel.: +43 316 873 7810. E-mail address: [email protected] (R. Prieler). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.08.056 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.
Radiative and convective heat transfer as well as temperature distributions in the gas phase and within the steel itself can be determined throughout the entire furnace domain, where conventional measurement techniques are of limited value only for such high temperature processes. For example, Jaklic et al. [2] investigated the transient cooling of steel billets during the transportation from the reheating furnace to the rolling mill and also considered the impact of the oxide scale layer on the heat transfer. A similar study on the impact of the oxide scale layer on the heat transfer was also done inside the furnace by Jang et al. [3]. Kim [4] predicted the radiative heat flux in a walking beam furnace with different absorption coefficients of the flue gas and emissivities of the slab. Although the convective heat flux was neglected, the slab temperature showed reasonably good agreement with the measurements. Due to the fact that wall and gas temperatures were assumed for different zones inside the furnace, no conclusion on the overall furnace efficiency (heat losses through flue gas, . . .) and pollutant emissions could be done. Similar simulations were performed by Emadi et al. [5], Tan et al. [6] and Jang et al. [3] considering the heat transfer to the furnace loads. A case study to optimise the reheating process in a walking beam furnace was carried out by Jaklic et al. [7] regarding the space between
676
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
the steel billets. The optimum space between the billets was found for three different billet sizes using a three-temperature model to calculate the heat transfer in the furnace. Nevertheless, such approaches are not able to predict the heat losses of the furnace and important information about the fluid flow (assumed convective heat fluxes) and the gas phase combustion like flame temperature and species concentrations is neglected. Simulation of the gas phase combustion focuses not only on the prediction of the radiative heat transfer inside the furnace. For example, Bhuiyan and Naser [8–10] carried out several simulations to investigate the heat fluxes for different oxy-fuel conditions also considering co-firing of biomass. A review of gas phase combustion of co-firing biomass was recently published in [11]. The gas phase combustion taking temperature and species concentrations into consideration was calculated by solving the governing transport equations as well as the transient heating of the steel simultaneously. This approach is computationally demanding due to the transient formulation of all the governing transport equations. For example, Han et al. [12] simulated the heating of a slab inside a bench scale reheating furnace with an in-house code using finite volume method (FVM). The eddy dissipation model proposed by Magnussen and Hjertager [13] was used for the combustion modelling which is computationally expensive due to the requirement of four additional transport equations for the considered species. Later Han et al. [14,15] started to use a non-premixed combustion model with presumed PDF (probability density function) to reduce the computational time. The transport of the furnace load to the next position was done with a user-defined function (UDF) in the commercial software package FLUENT. A more detailed analysis about furnace efficiency is possible thanks to the additional information about the fluid flow. A similar study was carried out by Gu et al. [16] on slab heating in a walking beam furnace simulating a small domain of a furnace only. A comparison of different methods to predict heat fluxes in reheating furnaces and transient heating of the furnace load was done by Singh and Talukdar [17] and Morgado et al. [18]. Singh and Talukdar investigated 4 models for a walking beam furnace. Three models simulated the radiative heat flux with different simplifications on the furnace and simulation models. The fourth model also considered the convective heat flux by solving the transport equations and showed an improved temperature prediction for the slabs. Morgado et al. used a detailed transient model, solving the governing transport equations with a presumed PDF to predict radiation and convection in the furnace. The second model subdivided the furnace in different zones calculating only the radiative heat transfer. It was found, that the second model was not able to predict the heat flux in great detail, but the predicted average volume temperatures of the slabs, however, showed a deviation of a mere 3%. In contrast to the works above, an iterative solution procedure for the steady-state gas phase combustion of the furnace and the transient heating of the steel billets was applied on a walking hearth furnace for reheating of steel billets. Both simulations were conducted with the software package FLUENT. First, the steady state solution of the combustion process in the furnace was used to calculate the total heat fluxes to the billets. A steady flamelet model with a detailed chemical kinetic mechanism was used in this study. Prieler et al. [19] revealed that a skeletal mechanism published by Peeters [20] can predict the flame for different combustion environments. Second, the transient heating characteristics of the billets was simulated using the calculated heat fluxes from the combustion model. The procedure was repeated until the temperature of the last billet did not change. This approach represents a novel method for coupling gas phase combustion and transient heating with future applicability for oxygen enhanced combustion due to the chemistry model used.
2. Furnace configuration To supply the rolling mill with steel billets several types of continuous reheating furnaces can be used. As shown in Section 1 many investigations were done on walking beam type furnaces were the slabs or billets are placed on stationary beams. The movement of the load is achieved by walking beams which raises and transports the load toward discharging. Operation in a walking hearth furnace, which is used in this study, is similar to a walking beam furnace. The steel billets are placed on stationary refractory blocks (hearths) and transported by the walking hearths. Compared to walking beams lower energy losses by cooling water for the transport system is necessary in a walking hearth furnace and also skid marks from the beams can be avoided. Nevertheless, the temperature distribution in a walking hearth furnace is not uniform because the bottom side is not exposed to the hot flue gas. A more detailed description of continuous industrial furnaces is given in [21]. A schematic overview of the heating process in the investigated walking hearth furnace is given in Fig. 1. Fig. 2 shows the CFD model of the walking hearth furnace without walls. The furnace model with walls is displayed in Fig. 3. The thicknesses of the furnace walls are given in Table 1. Because the furnace is symmetric in terms of the furnace length, one half only is displayed in Fig. 2. The dimensions of the combustion chamber are 17 m 13 m 1.6 m with an average load of 64 steel billets. The furnace is subdivided into three zones each with a different number and configuration of the burners. Flat flame burners were used for the combustion process with an overall fuel input of 18.2 MW. The pre-heating zone consists of 8 burners and is separated from heating and soaking zone by a water-cooled wall, which passes through the entire width of the furnace. The heating and soaking zones comprise 24 and 16 burners respectively. Natural-gas was used as fuel for the combustion process. The air (oxidizer) was preheated up to 310 °C. The furnace was operated under overall fuel lean conditions with an equivalence ratio of 0.97. In Eq. (1) the equivalence ratio is defined, where mfuel is the mass of fuel, mox is the mass of oxidizer and the subscript stands for stoichiometric conditions. The operating conditions, which were determined in several test runs, are summarized in Table 2.
mfuel =mox U¼ mfuel =mox st
ð1Þ
The billets were made of low-alloy steel with a length of 12 m and a cross section of 0.12 m 0.12 m. A billet is supplied to the furnace with a temperature of 20 °C and a density of 7800 kg/m3. The thermal conductivity and specific heat of the steel are displayed in Fig. 4. These properties were calculated with JMatPro. The residence time of a billet at each position inside the furnace was 100 s. This leads to an overall residence time of 6400 s and a production rate of approximately 46,400 kg/h. 3. Numerical models 3.1. Solution procedure The aim of this study was the coupling of the gas-phase combustion and the unsteady heating characteristic of the billets by CFD. The transport phenomena are periodically transient due to the movement of the billets inside the furnace as well as the charging and discharging process. An unsteady simulation of the gas phase and the heat transfer in the billets must thus be performed. For example Han et al. [14] performed such an analysis for a walking beam type furnace with approximately 560,000 cells. In many industrial furnaces, complicated burner geometry leads to fine
677
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
Fig. 1. Schematic figure of the furnace operation.
Fig. 2. Configuration of the walking hearth furnace without walls.
Table 1 Wall thickness of the furnace. Top wall (m) Side wall (m) Bottom wall (m)
0.29 0.47 0.6
Table 2 Operating conditions of the walking hearth furnace. Temperature air (°C) Temperature NG (°C) NG pre-heating zone (kg/h) NG heating zone (kg/h) NG soaking zone (kg/h) Equivalence ratio O2 in the flue gas (vol%) Fig. 3. Configuration of the walking hearth furnace with walls (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
grids with high a number of cells. Hence, an unsteady simulation is not applicable in many cases due to the high calculation time. An iterative solution procedure was thus chosen for combustion modelling and transient heat transfer in the billets for the work in this paper. In Fig. 5 the scheme of the iterative solution approach is displayed. The procedure was carried out with two alternating kinds of simulations:
310 25 290 780 250 0.97 0.7
Steady-state calculation of the gas phase combustion (see Section 3.2). Transient simulation of the heat transfer inside the billets (see Section 3.3). The procedure was started with a steady-state simulation of the gas phase combustion in the furnace. A temperature boundary condition for the billets in the furnace was assumed by a linear trend with respect to the furnace length (see Fig. 6). This was the initial condition for the procedure. Heat fluxes to the billets were
678
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
the arithmetic average of the linear temperature profile and the results of the transient simulations. The average temperature was used as the new boundary condition for the furnace simulation. Calculated temperatures of the billet were compared to the previous iteration step after each transient simulation. The relative temperature deviation between the most recent iteration step and the previous one was calculated using Eq. (2), where E is the relative deviation, ti;rec and t i;prev are the results from the most recent simulation and the previous iteration step (temperature in K). Convergence was determined by an average relative deviation E for all cells in the transient simulation below 0.5% which represents the standard tolerance of a type B thermocouple in this temperature range. Fig. 4. Thermal conductivity and specific heat of the steel depending on the temperature.
extracted from this furnace simulation and used as boundary condition for the transient simulation of the heat transfer in the billets. This simulation computed new surface temperatures ti for the billets depending on the position in the furnace. Following on from this the next iteration step was started by the definition of a new temperature boundary condition for the furnace simulation. This temperature profile was calculated by
E¼
jt i;rec t i;prev j ti;prev
ð2Þ
Furthermore, a solution procedure was tested without average temperature profiles. Temperature profiles from the transient calculation were imported directly to the next step. In this case, unreliable high temperatures (above melting point) were detected at the edges and corners of the billet after the second iteration step. This is a consequence of the strong dependence of the material properties on the temperature (see Fig. 4). Hence, the average temperature was introduced to the process and can be seen as a sort of under relaxation with a value of 0.5.
Fig. 5. Scheme of the solution procedure.
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
679
Fig. 6. Assumed temperature profile at the billet surfaces for the solution procedure.
3.2. Furnace simulation 3.2.1. Flow modelling The fluid flow in the furnace was calculated using the Reynolds-averaged Navier–Stokes equations given in Eqs. (3) and (4). For this purpose a double precision pressure-based solver was used for the steady-state solution of the continuity and momentum equation. For the simulations, the flue gas was assumed to be incompressible. Changes in the density due to the temperature were calculated by the ideal gas law. Operating pressure in the CFD model was set to 1013.25 mbar. SIMPLE and PRESTO! algorithms were applied for pressure–velocity coupling and pressure interpolation, respectively.
@ ðqui Þ ¼ 0 @xi @ @p @ @u @u 2 @u ðqui uj Þ ¼ þ l i þ j dij l @xj @xi @xj @xj @xi 3 @xl @ þ qu0i u0j @xj
ð3Þ
ð4Þ
For turbulence modelling a number of models are available in commercial CFD codes optimised for different applications. For example Bhuiyan et al. [22] tested the standard k-epsilon, RNG k-epsilon and k-omega models for a wavy fin-and-tube heat exchanger. A comparison with the experiments showed best agreement with the k-omega model. Based on this the k-omega model was used for case studies to investigate the effect of geometry of heat exchangers (see Bhuiyan et al. [23,24]). Nevertheless, the Reynolds stresses in the Eq. (4) were calculated in this study using two different turbulence models to test the impact on the solution especially in the vicinity of the burners (see Section 4.1) where high streamline curvature occurs. The turbulence models were tested with boundary conditions for the first iteration step, which means that a linear temperature profile at the billet surface was assumed. First, the realizable k-epsilon model proposed by Shih et al. [25] with a second order upwind scheme was tested. It had already been tested for several applications with axisymmetric jets and streamline curvatures (e.g. [26]). The same model was also investigated by Prieler et al. [19] for a jet flame in a lab-scale furnace and compared to the standard k-epsilon model and the Reynolds stress model (RSM). It was found, that realizable k-epsilon model predicted similar results as the RSM. Turbulent kinetic energy k and dissipation rate e were solved by just two additional transport equations (see [27]). Further information about the
transport equations and used model constants can be found in [25,27,28]. The second turbulence model was the Reynolds stress model (RSM) [29] combined with the QUICK scheme for discretization which leads to a better prediction of the anisotropic behaviour such as swirls and rapid changes on the strain rates and should be preferred for complex flows. In this case six additional transport equations for the Reynolds stresses had to be solved and increased the calculation time. For the first iteration step of the furnace simulation, a number of about 12,000 iterations was necessary to reach a converged solution. For the following furnace simulations 2000–3000 iterations only were sufficient due to the fact, that the solution of the previous iteration step was used as the initial condition for the simulation of the furnace. The convergence was determined by residuals for all equations below a value of 103. Temperature and species concentrations at the outlet as well as heat fluxes on the walls and billets were monitored for additional information about the simulation progress. Also the maximum error of the energy balance for all iteration steps of the furnace was 0.17% compared to the thermal input of 18.2 MW. 3.2.2. Boundary conditions and numerical grid Methane was used as fuel in the CFD model instead of naturalgas and air, with a composition of 21 vol% O2 and 79 vol% N2, as oxidizer. Whereas methane was supplied with 25 °C the oxidizer was pre-heated up to 310 °C. All inlets at the burners were modelled by mass-flow-inlet conditions with a turbulent intensity of 10%. Hydraulic diameters were set to 0.11 m and 0.0155 m for the fuel and oxidizer inlets. The mass-flow-rates for each zone as well as the equivalence ratio during furnace operation are given in Table 2. An oxygen volume fraction of 0.7 vol% at the outlet was detected due to fuel-lean conditions. A convective boundary condition was applied at the outer walls with a free stream temperature of 30 °C and a heat transfer coefficient of 20 W/(m2 K). The temperature boundary condition for the first iteration step was assumed by a linear profile according to Fig. 6. After the first iteration step, the temperature boundary condition was calculated from the results of the transient simulation of the billets. The emissivity of the billets and the furnace walls was defined with 0.5 and 0.7 respectively. The emissivity of the billets and the walls were determined based on values used in [6,15,18]. These authors used emissivities of 0.5 and 0.6 for the load and values between 0.6 and 0.75 for the walls. Because the furnace was symmetric along the length, only one half was modelled to reduce the number of cells. The
680
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
computational grid for the furnace simulation was made of an overall number of cells of 5,797,749 (tetrahedrons and hexahedrons) including combustion chamber and furnace walls (see Fig. 3) with a maximum edge length of approximately 400 mm in the region between flue gas outlet and pre-heating zone. These relatively large cells can be used in this region because low gradients of temperature, velocity and species concentrations were expected there. Furthermore, a simulation on a mesh with 8,026,188 cells was performed as a reference. The comparison of the results showed negligible differences on the temperature and species concentrations especially in the near burner region. Also the predicted heat fluxes were identical. Therefore, the smaller grid was used for further simulations. In Fig. 7 the combustion chamber of the smaller grid is shown with 4,830,514 cells only (without the solid walls). About 1.9 million cells were used for modelling the gas phase inside of the 24 flat flame burners to reach a high accuracy on the flame shape, temperature and species distribution inside the burners. The minimum edge length was 2.4 mm in the burner where the streamline curvatures were expected. In the conical part of the burner the edge length of the cells was between 9.3 and 32.6 mm compared to a burner diameter between 244 and 600 mm in this region. The solid walls consisted of approx. 1 million cells.
The turbulent flame is assumed to be an ensemble of small laminar diffusion flames. In such counter-flow diffusion flames the thermochemical state at each point between the fuel and oxidizer inlet was calculated before the CFD simulation. The one-dimensional analysis of the counter-flow diffusion flame was done with a skeletal reaction mechanism proposed by Peeters [20] using 17 species and 25 reversible reactions. The equations for the flamelet calculation can be found in [27]. Main advantage of this approach is that the chemical kinetic was considered without solving the complex flow field. The settings for the flamelet calculations are given in Table 3. The transport equations mentioned in Section 3.2.1 were solved for the Favre-averaged values of the fluctuating scalars, e.g. temperature and species concentrations. Thus, the relationship of instantaneous and average values of the scalars had to be considered by a probability density function (PDF) which describes the turbulence-chemistry interaction. A presumed b-probability density function p(f) was used for this investigation which can be thought as the fraction of time the fluid is in the state f [27]. The average value of the scalar /i , e.g. temperature, can be calculated by Eq. (6), where H is the mean enthalpy and /i is the instantaneous value of the scalar. It can be seen that the average value is dependent on the mean mixture fraction and the mean enthalpy. 02
3.2.3. Combustion and radiation models The aim of this study was to simulate the combustion process and the transient heating of the steel billets with high accuracy and a low calculation time. Hence the selection of the appropriate combustion and radiation models was crucial. Prieler et al. [19,28] investigated the eddy dissipation model, eddy dissipation concept model and steady flamelet model for different oxygen concentrations in the oxidizer. A minimum of the calculation time was found using the steady laminar flamelet model (SFM), although a detailed reaction mechanism was used in these studies. The chemistry calculations for the detailed mechanism were done before the CFD simulation and computationally demanding integrations for the chemical kinetic could be avoided during the simulation. The thermochemical state (mass fractions, density and temperature) in a laminar counter-flow diffusion flame is related only to the mean mixture fraction f [30,31]. The mean mixture fraction is defined by Eq. (5), where Z i , Z i;ox and Z i;fuel are the elemental mass fractions of element i in oxidizer and fuel.
f ¼
Z i Z i;ox Z i;fuel Z i;ox
ð5Þ
Also the mixture fraction variance f was used for integration of p(f) according to the presumed b-shape of the PDF [32]. The mixture fraction variance is defined in Eq. (7), where f and f are the instantaneous and the averaged mean mixture fraction. The calculation of the scalars, depending on the mixture fraction, mixture fraction variance and enthalpy, was done before the CFD simulation and the results were stored in look-up tables. Hence this approach reduced the chemistry to a simple mixing problem and kept the calculation time to a minimum, although a detailed reaction mechanism was used. All the chemistry calculations were preprocessed without solving the complex flow field in the furnace.
Z /i ¼
1
pðf Þ/i ðf ; HÞdf
ð6Þ
0 0
f ¼f f
ð7Þ
As described the averaged scalar /i is related to the mixture fraction, mixture fraction variance and the enthalpy. Thus two additional transport equations for the mean mixture fraction and the mixture fraction variance (Favre-averaged values) only needed to be solved (see Eqs. (8) and (9)) in the CFD code. The transport equation for the mean enthalpy is given in Eq. (10). Based on the
Fig. 7. Numerical grid of the furnace (a) and the burner (b).
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688 Table 3 Settings for flamelet and PDF calculation.
4. Results and discussion 4.1. Turbulence models
Flamelet calculation Maximum number of flamelets Number of grid points in flamelet Initial scalar dissipation (1/s) Scalar dissipation step (1/s)
20 32 0.01 5
PDF calculation Number of mean mixture fraction points Number of mixture fraction variance points Minimum temperature (°C) Number of mean enthalpy points
40 40 20 41
results for the mixture fraction, mixture fraction variance and enthalpy, the scalars (temperature, . . .) can be found in the lookup table without integrating the chemistry during the simulation process. It was possible to use a skeletal mechanism with 17 species, including OH, H and O radicals, and 25 reversible reactions within a short calculation time. Additional transport equations for the species could be avoided. The settings used for the flamelet and PDF calculations are given in Table 3.
r q~ v f ¼ r
r q~ v f 02 ¼ r
r q~ vH ¼ r
lt rf rt
ð8Þ
2 lt 02 02 þ C g lt rf C d q f rf rt k
kt r H þ Sh cp
ð9Þ
ð10Þ
The radiative transport equation (RTE) was solved in this study by the discrete ordinates (DO) model. Model coefficients proposed by Smith et al. [33] were applied. This model considers the RTE for a finite number of discrete angles. The RTE was solved for 128 directions which was a discretization of 4 4 directions per octant. Further information on the calculation of the radiative heat transfer is given [27]. 3.3. Transient simulation of billet heating The transient simulation of the heating characteristic of the billet was done on a computational grid with 86,400 hexahedrons for the solid. Thermal conductivity and specific heat of the billet were implemented by user-defined functions (see Fig. 4) in the CFD and the density was 7800 kg/m3. The initial temperature of the billet in the furnace was 20 °C. Boundary conditions at the walls of the billet were defined by profiles of the heat fluxes created as explained in Section 3.1. Heat fluxes were updated after 100 s because they were moved to the next position in the furnace. User-defined functions were created for updating the heat fluxes. The energy equation only was considered during the transient analysis with a second order spatial discretization and second order implicit for the transient formulation. A fixed time stepping method was chosen with a time step size of 1 s with a total number of 6400 steps. Convergence was achieved for one time step when the residual of the energy equation was below a value of 1011. The calculation time for the billet was approximately 2 h. In Eq. (11) the energy equation for the steel billets is given, where q is the density, k is the thermal conductivity and T is the temperature.
@ ðqhÞ ¼ r ðkrTÞ @t
681
ð11Þ
Due to the periodically transient simulation of the gas phase combustion and the billet heating characteristic, the computational time was kept to a minimum by using the steady flamelet model although a detailed chemistry was applied. In contrast to jet flames, which have been investigated by many researchers including Han et al. [14], Morgado et al. [18] and Singh and Talukdar [17], flat flame burners with high streamline curvature were used in the present work. Hence, a detailed numerical grid had to be used to calculate temperature and species concentrations in the flame. The numerical grid of the burners consisted of 1.9 million cells. Emadi et al. [5] carried out their research on a similar furnace with flat flame burners by a zone model to calculate the radiative heat flux in the furnace. The convective heat flux was determined by an empirical correlation. In this study, radiation and convection were calculated by CFD. For high accuracy the turbulence modelling is a crucial part and therefore the realizable k-epsilon model as well as the Reynolds stress model (RSM) were tested for a burner in the heating zone. Fig. 8 shows the calculated temperatures in the vicinity of the burner. The RSM predicted a similar flame shape and temperature. In both cases the flame temperature rose above 1600 °C at the beginning of the conical part of the burner. CFD models showed that the flame is in contact with the burner wall as well as the ceiling which indicated a good swirl. The flame length in this study was indicated by the main reaction zone, where the highest radical concentrations occur (for example OH radicals in Fig. 10). A slightly longer flame was predicted by the RSM which means that the flame is smaller than the flame calculated by the realizable k-epsilon model. This is a consequence of the better prediction of the swirl with the RSM and is highlighted by the radical concentration in the flame. In Fig. 9 the calculated wall temperatures in the burner for the RSM (a) and realizable k-epsilon model (b) are displayed. The realizable k-epsilon model predicted higher temperatures of approx. 1300 °C at the conical part of the burner wall. Due to the longer flames in the simulation with the RSM the wall temperature is 30–50 K lower than with the realizable k-epsilon model. A slightly asymmetric temperature distribution was detected with the realizable k-epsilon model due to the flue gas stream which carries away the hot gas from the flame. The realizable k-epsilon predicted lower swirl intensity and a shorter flame than the RSM and therefore the temperature distribution can be disturbed easily by the flue gas. This leads to a more asymmetric temperature distribution with the realizable k-epsilon model. During furnace operation the wall temperature at the burner was measured with a quotient pyrometer (Keller Optix Q PT70) at different positions but all with the same radial distance from burner axis (see dashed line in Fig. 9). Due to the turbulent fluctuations the wall temperature was measured in a range between 1250 and 1360 °C which is in good accordance with the predicted values of both turbulence models. Therefore both models are applicable for the investigated burners and furnace. In Fig. 10 the OH mole fraction is displayed which illustrates the main reaction zone. Higher radical concentrations in the burner were detected with the realizable k-epsilon model due to the shorter flame. Therefore, higher temperatures were predicted in the burner and at the burner walls which leads to an increasing OH concentration. Also the heat fluxes were compared due to their impact on the transient simulations of the billets. The relative errors of the heat fluxes between the RSM and the realizable k-epsilon model are displayed in Fig. 11. It is defined by Eq. (12), where Eq is the relative
682
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
Fig. 8. Temperature plots with the RSM (a) and realizable k-epsilon model (b).
Fig. 9. Wall temperatures at the burner with RSM (a) and realizable k-epsilon model (b).
Fig. 10. Plots of the OH mole fraction in the burner with the RSM (a) and realizable k-epsilon model (b).
error of the heat flux, Q_ RSM is the heat flux calculated with RSM and Q_ RKE denotes the heat flux predicted with the realizable k-epsilon model. Highest deviation can be seen for the heat losses through the wall with an error of 0.82%. This is a total difference of approximately 35 kW between the turbulence models. The two tested turbulence models showed basically similar results in the context of the flame shape and temperature as well as the predicted heat fluxes. The realizable k-epsilon model was thus used for all further simulations in this paper because of the lower computational demand by solving just two additional equations.
Q_
RSM Q_ RKE
Eq ¼
Q_ RSM
ð12Þ
4.2. Convergence criterion for the solution procedure and calculation time The calculated temperatures and relative deviations for each iteration step are displayed in Figs. 12 and 13. In these figures one half only of the last billet in the furnace is shown (before discharge). The first transient simulation showed very high temperatures at the end of the billet and the symmetry plane. A maximum value of 1365 °C was detected at the billet edge because of the higher thermal resistivity. After the third iteration step the surface temperature of the billet becomes more homogeneous. The difference between the maximum and minimum temperature decreased to 81 K and with the fifth iteration step the solution procedure was
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
683
First furnace simulation showed the highest calculation time at about 4 to 5 days. Because this result was used as the initial condition for the next furnace simulation, the calculation time could be reduced to 1 day. One transient calculation of the billet needed 2 h and had a minor impact on the duration of the overall process. It depicts a major decrease on computational time although a detailed chemical reaction mechanism was used for a good prediction of the flame and the heat transfer to the billets. 4.3. Furnace simulation
Fig. 11. Relative errors of the heat fluxes between realizable k-epsilon model and RSM.
aborted. At the top surface the fifth iteration step predicted an average temperature of 1180 °C with a maximum of 1207 °C. The relative deviation between the first transient simulation and the second iteration step were higher than 5% over the entire surface area except for some small regions (see Fig. 13). Despite an assumed temperature profile being used for the first furnace simulation, convergence for the solution procedure was already reached after the fifth iteration step. The average relative deviation was calculated with a value of 0.39% (4.7 K) between the fourth and fifth iteration step. Maximum deviation was found to be 0.87% resulting in a temperature deviation of approximately 10 K which was detected at the upper surface between the burners. Highest heat fluxes occur at these positions because of higher velocities and convective heat flux. Besides the relative deviations for all cells, the average temperatures of the surface and the entire volume change with 0.1 and 0.2 K between the last two iteration steps. The aim of this work was to find a solution procedure to calculate the gas phase combustion and heating characteristic with a high accuracy and a minimum on calculation time. Han et al. [14], for example, performed a transient simulation for combustion and heating of steel slabs simultaneously on a computational grid with approx. 570,000 cells which was finished after 6 days (140 h). In the investigated furnace the flat flame burners alone had to be meshed with about 1.9 million cells to compute the flame temperature and species concentrations. The overall cell number was 5,797,749 for the entire model. Thus a coupled transient simulation according to Han et al. would exceed the computational demand. The presented solution procedure in this paper was carried out within 10 days on an Intel i7 CPU with 3.2 GHz and 32 GB RAM.
In this section the results of the fifth iteration step are presented. In Fig. 14 the predicted heat fluxes are displayed. The thermal input of natural-gas was 18.2 MW according to the operating conditions. Due to the pre-heating of the air up to 310 °C an additional input of 1.9 MW was determined. CFD simulation calculated a heat transfer rate to the billets of 10.48 MW which is 52.1% of the total thermal input. This value can be seen as the furnace efficiency. Heat losses of the flue gas and the walls were 5.37 MW and 4.24 MW respectively. But it must be mentioned, that the predicted heat losses at the walls also contains losses due to billet charging and discharging. Fig. 15 shows the temperature and velocity in the furnace along the four burner rows. In the pre-heating zone a gas temperature between 600 and 800 °C was detected and therefore the radiative heat flux to the billets is low in the pre-heating zone with a value between 15 and 39 kW/m2 (see Fig. 16). In the pre-heating zone the convective heat flux reached a value of 3 kW/m2. The maximum convective heat flux in the furnace to the billet occurs between pre-heating and heating zone. A reduction of the cross section due to the water-cooled wall caused an increase of the flue gas velocity which is illustrated in the contour plot of Fig. 15. Higher velocities in this region (between pre-heating and heating zone) constitute a convective heat flux of 5 kW/m2. A comparison with another furnace investigated by Morgado et al. [18] showed that radiative heat fluxes in the pre-heating zone can reach approximately 85 kW/m2 when gas temperature is above 1000 °C in the pre-heating zone. Emadi et al. [5] calculated even higher radiative heat fluxes of about 110 kW/m2 in a walking hearth type furnace with a zone model for the radiative heat transfer. After the billet entered the heating zone the radiative heat flux increased up to 57.5 kW/m2 which is about 95% of the total heat flux in this zone, because the gas temperature was in a range of 1200–1300 °C. Up to the end of the furnace the radiative heat flux reached a maximum of 98% of the total heat flux to the billet. Radiative heat flux in the entire furnace was determined with 93% compared the overall heat transfer at the billet surface. Hence,
Fig. 12. Surface temperature of the last billet for each iteration step.
684
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
Fig. 13. Calculated relative deviations between the iteration steps: (a) Basic vs. first iteration; (b) First vs. second iteration; (c) Second vs. third iteration; (d) Third vs. fourth iteration.
1.6 m). In Fig. 17 the temperature trend calculated by CFD is interrupted at 8.5 m due to the water-cooled wall between pre-heating and heating zone which passes through the furnace at the measurement height of 0.8 m. It was found that a good accordance between measurement and CFD was reached for the entire furnace length. A maximum deviation was detected at the water-cooled wall with 83 K. Basically, the calculated temperatures were too low for all measurement points because of the higher heat losses through the walls including the losses of charging and discharging of the billets. Nevertheless, the mean deviation to the measured temperatures was 42 K. Based on this study it was determined that the CFD model is able to predict the gas phase combustion with high accuracy to the measured values and low computational time due to the use of the steady flamelet approach. Fig. 14. Thermal inputs and heat fluxes in the furnace.
4.4. Heating characteristic of the steel billets it can be concluded that convection had a minor effect on the heating of the billets especially in the heating and soaking zones. In addition to the analysis of the heat fluxes in the furnace, a comparison of the measured and simulated wall temperatures in the furnace was carried out. Measurement was done with the quotient pyrometer for both side walls. The values illustrated in Fig. 17 are the average temperatures of both side walls. The temperature measurement on the wall was done at a height of 0.8 m (overall
Knowledge about the temperature difference inside the billets during the heating is crucial for further processing on the rolling mill. When the billet enters the rolling mill a homogeneous temperature distribution as well as the desired target temperature is important during manufacturing to maintain optimum product quality. Furthermore the maximum temperature difference should be low to avoid deformations of the billets in the furnace. In this
Fig. 15. Contour plots of the furnace: (a) Temperature and (b) Velocity.
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
685
Fig. 16. Convective and radiative heat fluxes on the billets depending on furnace length.
Fig. 17. Measured and calculated wall temperatures in the furnace.
section the temperature distribution in the billet was investigated at different positions in the furnace. The temperature trend was observed along five planes along the furnace. These planes are displayed in Fig. 18 named ‘‘Level A” to ‘‘Level E”. The planes B and D (blue) were placed along the second and third row of burners whereas the other planes (grey) were between the burners. Nine temperature monitors labelled ‘‘t1” to ‘‘t9” were arranged inside the billet for each plane in the CFD model. Also a temperature monitor was set in the middle of the top surface (‘‘t_surf”). The calculated results for ‘‘Level A” only are displayed in Fig. 19, because the nine temperature monitors in all 5 planes showed the same trend along the furnace length. It was found that in the first half of the pre-heating zone the temperature inside the billet is quite homogenous with a small temperature difference. The maximum temperature difference in the billet until the end of the preheating zone increased to 99.9 K between t1 and t6. During the entire process t1 and t6 showed the highest temperature gap. The highest heat fluxes were detected on entering the heating zone and this leads to a faster increase in the temperature, especially for
the top monitors such as t1, t4 and t7. The lowest temperature occurs at position t6 where the heat transfer took place mainly due to conduction from the furnace bottom. The temperature difference reached a global maximum at 12.7 m. At this position the temperature at t6 started to increase rapidly. This is a result of the drop on the specific heat capacity (see Fig. 4) when the temperature is above 830 °C. For this condition the specific heat dropped from over 1000 to 600 J/(kgK) and the position t6 heats up faster than t1. In Fig. 20 the trend of the temperature difference between t1 and t6 is illustrated at all observed planes. Inside the billets the temperature difference started to increase with a constant rate in all planes. At approximately 4 m, when the billets reached the first burners of the pre-heating zone, the difference at the levels A and B were about 15 K lower than the other planes causing a slightly higher heat flux in the middle of the billets. As mentioned above, the temperature difference was determined at approximately 100 K when entering the heating zone. Due to the vast increase in the radiative heat flux, the temperature difference inside the steel increased faster until a maximum of 192 K was reached at
686
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
Fig. 18. Monitor planes through the furnace (grey: between the burners; blue: under the burners). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 19. Temperature trend in the billet at Level A.
Fig. 20. Temperature difference between t1 and t6 for all monitor planes (Level A to E).
level C, whereas a value of 170 K was observed at level B. A first decrease on the temperature difference in the billets was at 13 m which is still part of the heating zone. From this point to the end
of the soaking zone, the temperature differences dropped to a value between 42 and 49 K. Material properties with high dependency on the temperature can be assumed as a crucial part for
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
687
Fig. 21. Billet surface temperature at different positions: (a) 4.9 m; (b) 7.6 m; (c) 10.3 m; (d) 13.8 m; (e) 15.1 m; (f) 17 m.
modelling the heating characteristic and predicting the temperature gradients in the billets. In addition to the temperatures inside the furnace, the temperature trend on the upper surface (labelled with ‘‘t_surf”) was also considered in this study. At the beginning of the pre-heating zone the highest temperature difference occurred at 2.4 m with a value of 25.5 K between the levels A and E due to a higher convective heat flux to the billet. Increasing radiative heat flux in the heating and soaking zone leads to a more homogeneous temperature distribution on the surface. As a consequence the temperature difference becomes lower than 15 K in the heating zone. At the end of the furnace a temperature difference at the billet surface was determined with a maximum value of 7.6 K between level B and C. Contour plots of the surface temperature at different locations inside the furnace are shown in Fig. 21. The position at 4.9 m is below the first burners in the pre-heating zone. At this location CFD predicted higher temperatures at the corners and edges of the billets due to the higher thermal resistivity. A homogeneous surface temperature was determined at the end of the preheating zone due to increasing radiative heat transfer. This is illustrated by Fig. 21(b), which shows the billet at the end of the preheating zone. A higher temperature gradient was detected only at the corners of the billet. The same phenomenon was found at the end of the pre-heating zone at 10.3 m. After the billet was supported to the heating zone the gradient decreases at the edges due to higher exposure on the radiation as it can be seen in Fig. 21(d) at 13.8 m. This billet was located below the first burner row in the heating zone. In the soaking zone a further improvement in the temperature distribution was determined by the CFD simulation. Fig. 21(e) and (f) illustrate the billet at the beginning and the end of the soaking zone where a further decrease on the temperature difference at the surface leads to high product quality. It was pointed out that the higher radiative heat flux had a crucial impact on temperature distribution on the surface, which should be optimised. This could be achieved in future by adapting the furnace for oxygen enriched combustion and an investigation of the heating characteristics of the billets by CFD simulation.
5. Conclusion In this study a novel iterative solution procedure was tested to consider the transient heating process of steel billets in a walking hearth furnace. The iteration was carried out in two steps with a steady-state simulation of the gas phase combustion and the transient heat transfer in the billets. This approach was found to be a good alternative to zone models which calculate the radiative heat flux and convection by assumed zone temperatures and empirical correlations for the heat transfer coefficient (e.g. Kim et al. [4]) or
computationally demanding simultaneous transient simulations of gas phase combustion and heating of the load (e.g. Han et al. [14]). The overall calculation time for the combustion and transient heating was 10 days with consideration of radiative and convective heat transfer. The results of the steady-state simulation predicted that 52.1% (10.48 MW) of the total thermal input was transported to the steel billets. The overall radiative heat transfer to the billets was 93% of these 10.48 MW. The measured wall temperatures inside the furnace showed good accordance to the predicted values with an average deviation of 42 K for the entire furnace length. A constant increase of the temperature difference inside the billets was detected in the pre-heating zone with a value of 99.9 K at the end. The maximum was reached inside the billets in the heating zone where a peak temperature difference of 192 K was reached. This can be explained by the highest heat transfer rate in the furnace which leads to an increased temperature gradient. From this peak value the gradient decreased until the end of the furnace between 42 and 49 K. Optimisation on the furnace operation can be done without expensive test runs. This investigation forms the basis for the future simulation of the furnace under oxygen enriched conditions up to 25 vol% O2 to increase furnace efficiency with low modifications on the furnace. Due to the low enrichment level the product quality is expected to be the same. For this account simulations will be carried out to determine that oxy-fuel combustion can be applied without damage to the furnace or detriment to the desired product quality (temperature distribution). Conflict of interest There is no conflict of interest. Acknowledgment This work was financially supported by the Austrian Research Promotion Agency (FFG), ‘‘Competence-Headquarter: Ökologische und ökonomische Optimierung von industriellen Hochtemperaturprozessen durch Sauerstoffeinsatz” (Project 848281, eCall 4944044). References [1] Intergovernmental Panel on Climate Change (IPPC). IPPC special report on carbon dioxide capture and storage, Cambridge University Press, Cambridge, UK and New York, NY, USA, 2005. [2] A. Jaklic, B. Glogovac, T. Kolenko, B. Zupancic, B. Tezak, A simulation of heat transfer during billet transport, Appl. Therm. Eng. 22 (2002) 873–883. [3] J.H. Jang, D.E. Lee, M.Y. Kim, H.G. Kim, Investigation of the slab heating characteristics in a reheating furnace with the formation and growth of scale on the slab surface, Int. J. Heat Mass Transfer 53 (2010) 4326–4332.
688
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
[4] M.Y. Kim, A heat transfer model for the analysis of transient heating of the slab in a direct-fired walking beam type reheating furnace, Int. J. Heat Mass Transfer 50 (2007) 3740–3748. [5] A. Emadi, A. Saboonchi, M. Taheri, S. Hassanpour, Heating characteristics of billet in a walking hearth type reheating furnace, Appl. Therm. Eng. 63 (2014) 396–405. [6] C.K. Tan, J. Jenkins, J. Ward, J. Broughton, A. Heeley, Zone modelling of the thermal performances of a large-scale bloom reheating furnace, Appl. Therm. Eng. 50 (2013) 1111–1118. [7] A. Jaklic, T. Kolenko, B. Zupancic, The influence of the space between the billets on the productivity of a continuous walking-beam furnace, Appl. Therm. Eng. 25 (2005) 783–795. [8] A.A. Bhuiyan, J. Naser, Numerical modeling of biomass co-combustion with pulverized coal in a small scale furnace, Procedia Eng. 105 (2015) 504–511. [9] A.A. Bhuiyan, J. Naser, Numerical modelling of oxy fuel combustion, the effect of radiative and convective heat transfer and burnout, Fuel 139 (2015) 268–284. [10] A.A. Bhuiyan, J. Naser, Computational modelling of co-firing of biomass with coal under oxy-fuel condition in a small scale furnace, Fuel 143 (2015) 455–466. [11] M.R. Karim, J. Naser, Progress in numerical modelling of packed bed biomass combustion, in: 19th Australasian Fluid Mechanics Conference, Melbourne, Australia, 2014. [12] S.H. Han, S.W. Baek, S.H. Kang, C.Y. Kim, Numerical analysis of heating characteristics of a slab in a bench scale reheating furnace, Int. J. Heat Mass Transfer 50 (2007) 2019–2023. [13] B.F. Magnussen, B.H. Hjertager, On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion, in: 16th Symposium (Int) on Combustion, The Combust Inst, 1976. [14] S.H. Han, D. Chang, C.Y. Kim, A numerical analysis of slab heating characteristics in a walking beam type reheating furnace, Int. J. Heat Mass Transfer 53 (2010) 3855–3861. [15] S.H. Han, D. Chang, Optimum residence time analysis for a walking beam type reheating furnace, Int. J. Heat Mass Transfer 55 (2012) 4079–4087. [16] M.Y. Gu, G. Chen, X. Liu, C. Wu, H. Chu, Numerical simulation of slab heating process in a regenerative walking beam reheating furnace, Int. J. Heat Mass Transfer 76 (2014) 405–410. [17] V.K. Singh, P. Talukdar, Comparisons of different heat transfer models of a walking beam type reheat furnace, Int. Commun. Heat Mass Transfer 47 (2013) 20–26.
[18] T. Morgado, P.J. Coelho, P. Talukdar, Assessment of uniform temperature assumption in zoning on the numerical simulation of a walking beam reheating furnace, Appl. Therm. Eng. 76 (2015) 496–508. [19] R. Prieler, M. Demuth, D. Spoljaric, C. Hochenauer, Numerical investigation of the steady flamelet approach under different combustion environments, Fuel 140 (2015) 731–743. [20] T.W.J. Peeters, Numerical modeling oft turbulence natural-gas diffusion flames (Ph.D. thesis), Delft Technical University, Delft, The Netherlands, 1995. [21] W. Trinks, Industrial Furnaces, sixth ed., Jon Wiley & Sons, Inc., Hoboken, New Jersey, 2004. [22] A.A. Bhuiyan, A.K.M. Sadrul Islam, M. Ruhul Amin, Numerical study of 3D thermal and hydraulic characteristics of wavy fin and tube heat exchanger, Front. Heat Mass Transfer 3 (2012). [23] A.A. Bhuiyan, M. Ruhul Amin, J. Naser, A.K.M. Sadrul Islam, Effects of geometric parameters for wavy finned-tube heat exchanger in turbulent flow: a CFD modeling, Front. Heat Mass Transfer 6 (2015). [24] A.A. Bhuiyan, M. Ruhul Amin, R. Karim, A.K.M. Sadrul Islam, Plate fin and tube heat exchanger modeling: effects of performance parameters for turbulent flow regime, Int. J. Automot. Mech. Eng. 9 (2014) 1768–1781. [25] T.-H. Shih, W.W. Liou, A. Shabbir, Z. Yang, J. Zhu, A new k-e eddy viscosity model for high reynolds number turbulent flows – model development and validation, Comput. Fluids 24 (3) (1995) 227–238. [26] S.-E. Kim, D. Choudhury, B. Patel, Computations of complex turbulent flows using the commercial code ANSYS FLUENT, in: Proceedings of the ICASE/LaRC/ AFOSR Symposium on Modeling Complex Turbulent Flows, Hampton, Virginia, 1997. [27] ANSYS Fluent Theory Guide 13.0, ANSYS, Inc, Canonsburg, PA, USA, 2010. [28] R. Prieler, M. Demuth, D. Spoljaric, C. Hochenauer, Evaluation of a steady flamelet approach for use in oxy-fuel combustion, Fuel 118 (2014) 55–68. [29] B.E. Launder, G.J. Reece, W. Rodi, Progress in the development of a Reynoldsstress turbulence closure, J. Fluid Mech. 68 (3) (1975) 537–566. [30] N. Peters, Laminar diffusion flamelet models in non-premixed turbulent combustion, Prog. Energy Combust. Sci. 10 (1984) 319–339. [31] N. Peters, Laminar flamelet concepts in turbulent combustion, in: 21st Symposium (Int) on Combustion, The Combust Inst, 1986, pp. 1231–1250. [32] W.P. Jones, J.H. Whitelaw, Calculation methods for reacting turbulent flows: a review, Combust. Flame 48 (1982) 1–26. [33] T.F. Smith, Z.F. Shen, J.N. Friedman, Evaluation of coefficients for the weighted sum of gray gases model, J. Heat Transfer 104 (1982) 602–608.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Prediction of the heating characteristic of billets in a walking hearth type reheating furnace using CFD Rene Prieler a,⇑, Bernhard Mayr a, Martin Demuth b, Burkhardt Holleis b, Christoph Hochenauer a a b
Institute of Thermal Engineering, Graz University of Technology, Inffeldgasse 25/B, A-8010 Graz, Austria Messer Austria GmbH – Kompetenzzentrum Metallurgie, Industriestraße 5, A-2352 Gumpoldskirchen, Austria
a r t i c l e
i n f o
Article history: Received 17 June 2015 Received in revised form 3 August 2015 Accepted 4 August 2015 Available online 26 September 2015 Keywords: Combustion modelling Reheating furnace Computational fluid dynamics Heating characteristic Steady flamelet
a b s t r a c t A natural-gas fired walking hearth type furnace for the reheating of steel billets was investigated in this work. A novel numerical approach was used to predict the gas phase combustion, heat transfer and transient heating characteristics of the billets in the furnace. In contrast to conventional coupled simulation of the heat transfer and the heating characteristics, the iterative solution procedure used considers steady-state gas phase combustion of the furnace and unsteady simulation of the billets separately. Due to the steady flamelet model for combustion modelling the calculation time was kept to a minimum although a detailed reaction mechanism with 17 species and 25 reversible reactions was used. The reaction mechanism is applicable for different combustion environments and can be used for future investigations of the furnace after adaption for oxygen enriched combustion. Burner geometry was modelled in detail for a high accuracy prediction of the flame shape of the flat flame burners as well as temperature and species concentrations in the vicinity of the burner. The simulation revealed that 93% of the total heat flux to the billets was contributed by radiation. Transient heating showed a maximum temperature difference inside the billets in the heating zone at 192 K decreasing to a value between 42 and 49 K at the end of the furnace. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Large scale reheating furnaces are employed in the steel industry for heating up steel billets or slabs to a desired temperature level for further treatment in a rolling mill. Furthermore, a uniform temperature distribution is desired inside the billets with optimised operating conditions regarding fuel consumption, residence time, NOx and CO2 emissions etc. According to an IPCC report, about 60% of global CO2 emissions are produced by fossil fuels fired processes [1]. In the world today environmental issues and cost efficient production are the subject of major attention in numerous industrial sectors. Significant efforts are being made to increase efficiency against this background. In the past few years as a result of increasing computer power, many researchers have resorted to numerical methods for investigating reheating furnaces as these are applied in the steel manufacturing process. Numerical models can offer a detailed view of the transport phenomena in the furnaces and can be a helpful tool to optimise the heating process, transportation of the furnace loads and other phenomena. ⇑ Corresponding author. Tel.: +43 316 873 7810. E-mail address: [email protected] (R. Prieler). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.08.056 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.
Radiative and convective heat transfer as well as temperature distributions in the gas phase and within the steel itself can be determined throughout the entire furnace domain, where conventional measurement techniques are of limited value only for such high temperature processes. For example, Jaklic et al. [2] investigated the transient cooling of steel billets during the transportation from the reheating furnace to the rolling mill and also considered the impact of the oxide scale layer on the heat transfer. A similar study on the impact of the oxide scale layer on the heat transfer was also done inside the furnace by Jang et al. [3]. Kim [4] predicted the radiative heat flux in a walking beam furnace with different absorption coefficients of the flue gas and emissivities of the slab. Although the convective heat flux was neglected, the slab temperature showed reasonably good agreement with the measurements. Due to the fact that wall and gas temperatures were assumed for different zones inside the furnace, no conclusion on the overall furnace efficiency (heat losses through flue gas, . . .) and pollutant emissions could be done. Similar simulations were performed by Emadi et al. [5], Tan et al. [6] and Jang et al. [3] considering the heat transfer to the furnace loads. A case study to optimise the reheating process in a walking beam furnace was carried out by Jaklic et al. [7] regarding the space between
676
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
the steel billets. The optimum space between the billets was found for three different billet sizes using a three-temperature model to calculate the heat transfer in the furnace. Nevertheless, such approaches are not able to predict the heat losses of the furnace and important information about the fluid flow (assumed convective heat fluxes) and the gas phase combustion like flame temperature and species concentrations is neglected. Simulation of the gas phase combustion focuses not only on the prediction of the radiative heat transfer inside the furnace. For example, Bhuiyan and Naser [8–10] carried out several simulations to investigate the heat fluxes for different oxy-fuel conditions also considering co-firing of biomass. A review of gas phase combustion of co-firing biomass was recently published in [11]. The gas phase combustion taking temperature and species concentrations into consideration was calculated by solving the governing transport equations as well as the transient heating of the steel simultaneously. This approach is computationally demanding due to the transient formulation of all the governing transport equations. For example, Han et al. [12] simulated the heating of a slab inside a bench scale reheating furnace with an in-house code using finite volume method (FVM). The eddy dissipation model proposed by Magnussen and Hjertager [13] was used for the combustion modelling which is computationally expensive due to the requirement of four additional transport equations for the considered species. Later Han et al. [14,15] started to use a non-premixed combustion model with presumed PDF (probability density function) to reduce the computational time. The transport of the furnace load to the next position was done with a user-defined function (UDF) in the commercial software package FLUENT. A more detailed analysis about furnace efficiency is possible thanks to the additional information about the fluid flow. A similar study was carried out by Gu et al. [16] on slab heating in a walking beam furnace simulating a small domain of a furnace only. A comparison of different methods to predict heat fluxes in reheating furnaces and transient heating of the furnace load was done by Singh and Talukdar [17] and Morgado et al. [18]. Singh and Talukdar investigated 4 models for a walking beam furnace. Three models simulated the radiative heat flux with different simplifications on the furnace and simulation models. The fourth model also considered the convective heat flux by solving the transport equations and showed an improved temperature prediction for the slabs. Morgado et al. used a detailed transient model, solving the governing transport equations with a presumed PDF to predict radiation and convection in the furnace. The second model subdivided the furnace in different zones calculating only the radiative heat transfer. It was found, that the second model was not able to predict the heat flux in great detail, but the predicted average volume temperatures of the slabs, however, showed a deviation of a mere 3%. In contrast to the works above, an iterative solution procedure for the steady-state gas phase combustion of the furnace and the transient heating of the steel billets was applied on a walking hearth furnace for reheating of steel billets. Both simulations were conducted with the software package FLUENT. First, the steady state solution of the combustion process in the furnace was used to calculate the total heat fluxes to the billets. A steady flamelet model with a detailed chemical kinetic mechanism was used in this study. Prieler et al. [19] revealed that a skeletal mechanism published by Peeters [20] can predict the flame for different combustion environments. Second, the transient heating characteristics of the billets was simulated using the calculated heat fluxes from the combustion model. The procedure was repeated until the temperature of the last billet did not change. This approach represents a novel method for coupling gas phase combustion and transient heating with future applicability for oxygen enhanced combustion due to the chemistry model used.
2. Furnace configuration To supply the rolling mill with steel billets several types of continuous reheating furnaces can be used. As shown in Section 1 many investigations were done on walking beam type furnaces were the slabs or billets are placed on stationary beams. The movement of the load is achieved by walking beams which raises and transports the load toward discharging. Operation in a walking hearth furnace, which is used in this study, is similar to a walking beam furnace. The steel billets are placed on stationary refractory blocks (hearths) and transported by the walking hearths. Compared to walking beams lower energy losses by cooling water for the transport system is necessary in a walking hearth furnace and also skid marks from the beams can be avoided. Nevertheless, the temperature distribution in a walking hearth furnace is not uniform because the bottom side is not exposed to the hot flue gas. A more detailed description of continuous industrial furnaces is given in [21]. A schematic overview of the heating process in the investigated walking hearth furnace is given in Fig. 1. Fig. 2 shows the CFD model of the walking hearth furnace without walls. The furnace model with walls is displayed in Fig. 3. The thicknesses of the furnace walls are given in Table 1. Because the furnace is symmetric in terms of the furnace length, one half only is displayed in Fig. 2. The dimensions of the combustion chamber are 17 m 13 m 1.6 m with an average load of 64 steel billets. The furnace is subdivided into three zones each with a different number and configuration of the burners. Flat flame burners were used for the combustion process with an overall fuel input of 18.2 MW. The pre-heating zone consists of 8 burners and is separated from heating and soaking zone by a water-cooled wall, which passes through the entire width of the furnace. The heating and soaking zones comprise 24 and 16 burners respectively. Natural-gas was used as fuel for the combustion process. The air (oxidizer) was preheated up to 310 °C. The furnace was operated under overall fuel lean conditions with an equivalence ratio of 0.97. In Eq. (1) the equivalence ratio is defined, where mfuel is the mass of fuel, mox is the mass of oxidizer and the subscript stands for stoichiometric conditions. The operating conditions, which were determined in several test runs, are summarized in Table 2.
mfuel =mox U¼ mfuel =mox st
ð1Þ
The billets were made of low-alloy steel with a length of 12 m and a cross section of 0.12 m 0.12 m. A billet is supplied to the furnace with a temperature of 20 °C and a density of 7800 kg/m3. The thermal conductivity and specific heat of the steel are displayed in Fig. 4. These properties were calculated with JMatPro. The residence time of a billet at each position inside the furnace was 100 s. This leads to an overall residence time of 6400 s and a production rate of approximately 46,400 kg/h. 3. Numerical models 3.1. Solution procedure The aim of this study was the coupling of the gas-phase combustion and the unsteady heating characteristic of the billets by CFD. The transport phenomena are periodically transient due to the movement of the billets inside the furnace as well as the charging and discharging process. An unsteady simulation of the gas phase and the heat transfer in the billets must thus be performed. For example Han et al. [14] performed such an analysis for a walking beam type furnace with approximately 560,000 cells. In many industrial furnaces, complicated burner geometry leads to fine
677
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
Fig. 1. Schematic figure of the furnace operation.
Fig. 2. Configuration of the walking hearth furnace without walls.
Table 1 Wall thickness of the furnace. Top wall (m) Side wall (m) Bottom wall (m)
0.29 0.47 0.6
Table 2 Operating conditions of the walking hearth furnace. Temperature air (°C) Temperature NG (°C) NG pre-heating zone (kg/h) NG heating zone (kg/h) NG soaking zone (kg/h) Equivalence ratio O2 in the flue gas (vol%) Fig. 3. Configuration of the walking hearth furnace with walls (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
grids with high a number of cells. Hence, an unsteady simulation is not applicable in many cases due to the high calculation time. An iterative solution procedure was thus chosen for combustion modelling and transient heat transfer in the billets for the work in this paper. In Fig. 5 the scheme of the iterative solution approach is displayed. The procedure was carried out with two alternating kinds of simulations:
310 25 290 780 250 0.97 0.7
Steady-state calculation of the gas phase combustion (see Section 3.2). Transient simulation of the heat transfer inside the billets (see Section 3.3). The procedure was started with a steady-state simulation of the gas phase combustion in the furnace. A temperature boundary condition for the billets in the furnace was assumed by a linear trend with respect to the furnace length (see Fig. 6). This was the initial condition for the procedure. Heat fluxes to the billets were
678
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
the arithmetic average of the linear temperature profile and the results of the transient simulations. The average temperature was used as the new boundary condition for the furnace simulation. Calculated temperatures of the billet were compared to the previous iteration step after each transient simulation. The relative temperature deviation between the most recent iteration step and the previous one was calculated using Eq. (2), where E is the relative deviation, ti;rec and t i;prev are the results from the most recent simulation and the previous iteration step (temperature in K). Convergence was determined by an average relative deviation E for all cells in the transient simulation below 0.5% which represents the standard tolerance of a type B thermocouple in this temperature range. Fig. 4. Thermal conductivity and specific heat of the steel depending on the temperature.
extracted from this furnace simulation and used as boundary condition for the transient simulation of the heat transfer in the billets. This simulation computed new surface temperatures ti for the billets depending on the position in the furnace. Following on from this the next iteration step was started by the definition of a new temperature boundary condition for the furnace simulation. This temperature profile was calculated by
E¼
jt i;rec t i;prev j ti;prev
ð2Þ
Furthermore, a solution procedure was tested without average temperature profiles. Temperature profiles from the transient calculation were imported directly to the next step. In this case, unreliable high temperatures (above melting point) were detected at the edges and corners of the billet after the second iteration step. This is a consequence of the strong dependence of the material properties on the temperature (see Fig. 4). Hence, the average temperature was introduced to the process and can be seen as a sort of under relaxation with a value of 0.5.
Fig. 5. Scheme of the solution procedure.
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
679
Fig. 6. Assumed temperature profile at the billet surfaces for the solution procedure.
3.2. Furnace simulation 3.2.1. Flow modelling The fluid flow in the furnace was calculated using the Reynolds-averaged Navier–Stokes equations given in Eqs. (3) and (4). For this purpose a double precision pressure-based solver was used for the steady-state solution of the continuity and momentum equation. For the simulations, the flue gas was assumed to be incompressible. Changes in the density due to the temperature were calculated by the ideal gas law. Operating pressure in the CFD model was set to 1013.25 mbar. SIMPLE and PRESTO! algorithms were applied for pressure–velocity coupling and pressure interpolation, respectively.
@ ðqui Þ ¼ 0 @xi @ @p @ @u @u 2 @u ðqui uj Þ ¼ þ l i þ j dij l @xj @xi @xj @xj @xi 3 @xl @ þ qu0i u0j @xj
ð3Þ
ð4Þ
For turbulence modelling a number of models are available in commercial CFD codes optimised for different applications. For example Bhuiyan et al. [22] tested the standard k-epsilon, RNG k-epsilon and k-omega models for a wavy fin-and-tube heat exchanger. A comparison with the experiments showed best agreement with the k-omega model. Based on this the k-omega model was used for case studies to investigate the effect of geometry of heat exchangers (see Bhuiyan et al. [23,24]). Nevertheless, the Reynolds stresses in the Eq. (4) were calculated in this study using two different turbulence models to test the impact on the solution especially in the vicinity of the burners (see Section 4.1) where high streamline curvature occurs. The turbulence models were tested with boundary conditions for the first iteration step, which means that a linear temperature profile at the billet surface was assumed. First, the realizable k-epsilon model proposed by Shih et al. [25] with a second order upwind scheme was tested. It had already been tested for several applications with axisymmetric jets and streamline curvatures (e.g. [26]). The same model was also investigated by Prieler et al. [19] for a jet flame in a lab-scale furnace and compared to the standard k-epsilon model and the Reynolds stress model (RSM). It was found, that realizable k-epsilon model predicted similar results as the RSM. Turbulent kinetic energy k and dissipation rate e were solved by just two additional transport equations (see [27]). Further information about the
transport equations and used model constants can be found in [25,27,28]. The second turbulence model was the Reynolds stress model (RSM) [29] combined with the QUICK scheme for discretization which leads to a better prediction of the anisotropic behaviour such as swirls and rapid changes on the strain rates and should be preferred for complex flows. In this case six additional transport equations for the Reynolds stresses had to be solved and increased the calculation time. For the first iteration step of the furnace simulation, a number of about 12,000 iterations was necessary to reach a converged solution. For the following furnace simulations 2000–3000 iterations only were sufficient due to the fact, that the solution of the previous iteration step was used as the initial condition for the simulation of the furnace. The convergence was determined by residuals for all equations below a value of 103. Temperature and species concentrations at the outlet as well as heat fluxes on the walls and billets were monitored for additional information about the simulation progress. Also the maximum error of the energy balance for all iteration steps of the furnace was 0.17% compared to the thermal input of 18.2 MW. 3.2.2. Boundary conditions and numerical grid Methane was used as fuel in the CFD model instead of naturalgas and air, with a composition of 21 vol% O2 and 79 vol% N2, as oxidizer. Whereas methane was supplied with 25 °C the oxidizer was pre-heated up to 310 °C. All inlets at the burners were modelled by mass-flow-inlet conditions with a turbulent intensity of 10%. Hydraulic diameters were set to 0.11 m and 0.0155 m for the fuel and oxidizer inlets. The mass-flow-rates for each zone as well as the equivalence ratio during furnace operation are given in Table 2. An oxygen volume fraction of 0.7 vol% at the outlet was detected due to fuel-lean conditions. A convective boundary condition was applied at the outer walls with a free stream temperature of 30 °C and a heat transfer coefficient of 20 W/(m2 K). The temperature boundary condition for the first iteration step was assumed by a linear profile according to Fig. 6. After the first iteration step, the temperature boundary condition was calculated from the results of the transient simulation of the billets. The emissivity of the billets and the furnace walls was defined with 0.5 and 0.7 respectively. The emissivity of the billets and the walls were determined based on values used in [6,15,18]. These authors used emissivities of 0.5 and 0.6 for the load and values between 0.6 and 0.75 for the walls. Because the furnace was symmetric along the length, only one half was modelled to reduce the number of cells. The
680
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
computational grid for the furnace simulation was made of an overall number of cells of 5,797,749 (tetrahedrons and hexahedrons) including combustion chamber and furnace walls (see Fig. 3) with a maximum edge length of approximately 400 mm in the region between flue gas outlet and pre-heating zone. These relatively large cells can be used in this region because low gradients of temperature, velocity and species concentrations were expected there. Furthermore, a simulation on a mesh with 8,026,188 cells was performed as a reference. The comparison of the results showed negligible differences on the temperature and species concentrations especially in the near burner region. Also the predicted heat fluxes were identical. Therefore, the smaller grid was used for further simulations. In Fig. 7 the combustion chamber of the smaller grid is shown with 4,830,514 cells only (without the solid walls). About 1.9 million cells were used for modelling the gas phase inside of the 24 flat flame burners to reach a high accuracy on the flame shape, temperature and species distribution inside the burners. The minimum edge length was 2.4 mm in the burner where the streamline curvatures were expected. In the conical part of the burner the edge length of the cells was between 9.3 and 32.6 mm compared to a burner diameter between 244 and 600 mm in this region. The solid walls consisted of approx. 1 million cells.
The turbulent flame is assumed to be an ensemble of small laminar diffusion flames. In such counter-flow diffusion flames the thermochemical state at each point between the fuel and oxidizer inlet was calculated before the CFD simulation. The one-dimensional analysis of the counter-flow diffusion flame was done with a skeletal reaction mechanism proposed by Peeters [20] using 17 species and 25 reversible reactions. The equations for the flamelet calculation can be found in [27]. Main advantage of this approach is that the chemical kinetic was considered without solving the complex flow field. The settings for the flamelet calculations are given in Table 3. The transport equations mentioned in Section 3.2.1 were solved for the Favre-averaged values of the fluctuating scalars, e.g. temperature and species concentrations. Thus, the relationship of instantaneous and average values of the scalars had to be considered by a probability density function (PDF) which describes the turbulence-chemistry interaction. A presumed b-probability density function p(f) was used for this investigation which can be thought as the fraction of time the fluid is in the state f [27]. The average value of the scalar /i , e.g. temperature, can be calculated by Eq. (6), where H is the mean enthalpy and /i is the instantaneous value of the scalar. It can be seen that the average value is dependent on the mean mixture fraction and the mean enthalpy. 02
3.2.3. Combustion and radiation models The aim of this study was to simulate the combustion process and the transient heating of the steel billets with high accuracy and a low calculation time. Hence the selection of the appropriate combustion and radiation models was crucial. Prieler et al. [19,28] investigated the eddy dissipation model, eddy dissipation concept model and steady flamelet model for different oxygen concentrations in the oxidizer. A minimum of the calculation time was found using the steady laminar flamelet model (SFM), although a detailed reaction mechanism was used in these studies. The chemistry calculations for the detailed mechanism were done before the CFD simulation and computationally demanding integrations for the chemical kinetic could be avoided during the simulation. The thermochemical state (mass fractions, density and temperature) in a laminar counter-flow diffusion flame is related only to the mean mixture fraction f [30,31]. The mean mixture fraction is defined by Eq. (5), where Z i , Z i;ox and Z i;fuel are the elemental mass fractions of element i in oxidizer and fuel.
f ¼
Z i Z i;ox Z i;fuel Z i;ox
ð5Þ
Also the mixture fraction variance f was used for integration of p(f) according to the presumed b-shape of the PDF [32]. The mixture fraction variance is defined in Eq. (7), where f and f are the instantaneous and the averaged mean mixture fraction. The calculation of the scalars, depending on the mixture fraction, mixture fraction variance and enthalpy, was done before the CFD simulation and the results were stored in look-up tables. Hence this approach reduced the chemistry to a simple mixing problem and kept the calculation time to a minimum, although a detailed reaction mechanism was used. All the chemistry calculations were preprocessed without solving the complex flow field in the furnace.
Z /i ¼
1
pðf Þ/i ðf ; HÞdf
ð6Þ
0 0
f ¼f f
ð7Þ
As described the averaged scalar /i is related to the mixture fraction, mixture fraction variance and the enthalpy. Thus two additional transport equations for the mean mixture fraction and the mixture fraction variance (Favre-averaged values) only needed to be solved (see Eqs. (8) and (9)) in the CFD code. The transport equation for the mean enthalpy is given in Eq. (10). Based on the
Fig. 7. Numerical grid of the furnace (a) and the burner (b).
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688 Table 3 Settings for flamelet and PDF calculation.
4. Results and discussion 4.1. Turbulence models
Flamelet calculation Maximum number of flamelets Number of grid points in flamelet Initial scalar dissipation (1/s) Scalar dissipation step (1/s)
20 32 0.01 5
PDF calculation Number of mean mixture fraction points Number of mixture fraction variance points Minimum temperature (°C) Number of mean enthalpy points
40 40 20 41
results for the mixture fraction, mixture fraction variance and enthalpy, the scalars (temperature, . . .) can be found in the lookup table without integrating the chemistry during the simulation process. It was possible to use a skeletal mechanism with 17 species, including OH, H and O radicals, and 25 reversible reactions within a short calculation time. Additional transport equations for the species could be avoided. The settings used for the flamelet and PDF calculations are given in Table 3.
r q~ v f ¼ r
r q~ v f 02 ¼ r
r q~ vH ¼ r
lt rf rt
ð8Þ
2 lt 02 02 þ C g lt rf C d q f rf rt k
kt r H þ Sh cp
ð9Þ
ð10Þ
The radiative transport equation (RTE) was solved in this study by the discrete ordinates (DO) model. Model coefficients proposed by Smith et al. [33] were applied. This model considers the RTE for a finite number of discrete angles. The RTE was solved for 128 directions which was a discretization of 4 4 directions per octant. Further information on the calculation of the radiative heat transfer is given [27]. 3.3. Transient simulation of billet heating The transient simulation of the heating characteristic of the billet was done on a computational grid with 86,400 hexahedrons for the solid. Thermal conductivity and specific heat of the billet were implemented by user-defined functions (see Fig. 4) in the CFD and the density was 7800 kg/m3. The initial temperature of the billet in the furnace was 20 °C. Boundary conditions at the walls of the billet were defined by profiles of the heat fluxes created as explained in Section 3.1. Heat fluxes were updated after 100 s because they were moved to the next position in the furnace. User-defined functions were created for updating the heat fluxes. The energy equation only was considered during the transient analysis with a second order spatial discretization and second order implicit for the transient formulation. A fixed time stepping method was chosen with a time step size of 1 s with a total number of 6400 steps. Convergence was achieved for one time step when the residual of the energy equation was below a value of 1011. The calculation time for the billet was approximately 2 h. In Eq. (11) the energy equation for the steel billets is given, where q is the density, k is the thermal conductivity and T is the temperature.
@ ðqhÞ ¼ r ðkrTÞ @t
681
ð11Þ
Due to the periodically transient simulation of the gas phase combustion and the billet heating characteristic, the computational time was kept to a minimum by using the steady flamelet model although a detailed chemistry was applied. In contrast to jet flames, which have been investigated by many researchers including Han et al. [14], Morgado et al. [18] and Singh and Talukdar [17], flat flame burners with high streamline curvature were used in the present work. Hence, a detailed numerical grid had to be used to calculate temperature and species concentrations in the flame. The numerical grid of the burners consisted of 1.9 million cells. Emadi et al. [5] carried out their research on a similar furnace with flat flame burners by a zone model to calculate the radiative heat flux in the furnace. The convective heat flux was determined by an empirical correlation. In this study, radiation and convection were calculated by CFD. For high accuracy the turbulence modelling is a crucial part and therefore the realizable k-epsilon model as well as the Reynolds stress model (RSM) were tested for a burner in the heating zone. Fig. 8 shows the calculated temperatures in the vicinity of the burner. The RSM predicted a similar flame shape and temperature. In both cases the flame temperature rose above 1600 °C at the beginning of the conical part of the burner. CFD models showed that the flame is in contact with the burner wall as well as the ceiling which indicated a good swirl. The flame length in this study was indicated by the main reaction zone, where the highest radical concentrations occur (for example OH radicals in Fig. 10). A slightly longer flame was predicted by the RSM which means that the flame is smaller than the flame calculated by the realizable k-epsilon model. This is a consequence of the better prediction of the swirl with the RSM and is highlighted by the radical concentration in the flame. In Fig. 9 the calculated wall temperatures in the burner for the RSM (a) and realizable k-epsilon model (b) are displayed. The realizable k-epsilon model predicted higher temperatures of approx. 1300 °C at the conical part of the burner wall. Due to the longer flames in the simulation with the RSM the wall temperature is 30–50 K lower than with the realizable k-epsilon model. A slightly asymmetric temperature distribution was detected with the realizable k-epsilon model due to the flue gas stream which carries away the hot gas from the flame. The realizable k-epsilon predicted lower swirl intensity and a shorter flame than the RSM and therefore the temperature distribution can be disturbed easily by the flue gas. This leads to a more asymmetric temperature distribution with the realizable k-epsilon model. During furnace operation the wall temperature at the burner was measured with a quotient pyrometer (Keller Optix Q PT70) at different positions but all with the same radial distance from burner axis (see dashed line in Fig. 9). Due to the turbulent fluctuations the wall temperature was measured in a range between 1250 and 1360 °C which is in good accordance with the predicted values of both turbulence models. Therefore both models are applicable for the investigated burners and furnace. In Fig. 10 the OH mole fraction is displayed which illustrates the main reaction zone. Higher radical concentrations in the burner were detected with the realizable k-epsilon model due to the shorter flame. Therefore, higher temperatures were predicted in the burner and at the burner walls which leads to an increasing OH concentration. Also the heat fluxes were compared due to their impact on the transient simulations of the billets. The relative errors of the heat fluxes between the RSM and the realizable k-epsilon model are displayed in Fig. 11. It is defined by Eq. (12), where Eq is the relative
682
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
Fig. 8. Temperature plots with the RSM (a) and realizable k-epsilon model (b).
Fig. 9. Wall temperatures at the burner with RSM (a) and realizable k-epsilon model (b).
Fig. 10. Plots of the OH mole fraction in the burner with the RSM (a) and realizable k-epsilon model (b).
error of the heat flux, Q_ RSM is the heat flux calculated with RSM and Q_ RKE denotes the heat flux predicted with the realizable k-epsilon model. Highest deviation can be seen for the heat losses through the wall with an error of 0.82%. This is a total difference of approximately 35 kW between the turbulence models. The two tested turbulence models showed basically similar results in the context of the flame shape and temperature as well as the predicted heat fluxes. The realizable k-epsilon model was thus used for all further simulations in this paper because of the lower computational demand by solving just two additional equations.
Q_
RSM Q_ RKE
Eq ¼
Q_ RSM
ð12Þ
4.2. Convergence criterion for the solution procedure and calculation time The calculated temperatures and relative deviations for each iteration step are displayed in Figs. 12 and 13. In these figures one half only of the last billet in the furnace is shown (before discharge). The first transient simulation showed very high temperatures at the end of the billet and the symmetry plane. A maximum value of 1365 °C was detected at the billet edge because of the higher thermal resistivity. After the third iteration step the surface temperature of the billet becomes more homogeneous. The difference between the maximum and minimum temperature decreased to 81 K and with the fifth iteration step the solution procedure was
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
683
First furnace simulation showed the highest calculation time at about 4 to 5 days. Because this result was used as the initial condition for the next furnace simulation, the calculation time could be reduced to 1 day. One transient calculation of the billet needed 2 h and had a minor impact on the duration of the overall process. It depicts a major decrease on computational time although a detailed chemical reaction mechanism was used for a good prediction of the flame and the heat transfer to the billets. 4.3. Furnace simulation
Fig. 11. Relative errors of the heat fluxes between realizable k-epsilon model and RSM.
aborted. At the top surface the fifth iteration step predicted an average temperature of 1180 °C with a maximum of 1207 °C. The relative deviation between the first transient simulation and the second iteration step were higher than 5% over the entire surface area except for some small regions (see Fig. 13). Despite an assumed temperature profile being used for the first furnace simulation, convergence for the solution procedure was already reached after the fifth iteration step. The average relative deviation was calculated with a value of 0.39% (4.7 K) between the fourth and fifth iteration step. Maximum deviation was found to be 0.87% resulting in a temperature deviation of approximately 10 K which was detected at the upper surface between the burners. Highest heat fluxes occur at these positions because of higher velocities and convective heat flux. Besides the relative deviations for all cells, the average temperatures of the surface and the entire volume change with 0.1 and 0.2 K between the last two iteration steps. The aim of this work was to find a solution procedure to calculate the gas phase combustion and heating characteristic with a high accuracy and a minimum on calculation time. Han et al. [14], for example, performed a transient simulation for combustion and heating of steel slabs simultaneously on a computational grid with approx. 570,000 cells which was finished after 6 days (140 h). In the investigated furnace the flat flame burners alone had to be meshed with about 1.9 million cells to compute the flame temperature and species concentrations. The overall cell number was 5,797,749 for the entire model. Thus a coupled transient simulation according to Han et al. would exceed the computational demand. The presented solution procedure in this paper was carried out within 10 days on an Intel i7 CPU with 3.2 GHz and 32 GB RAM.
In this section the results of the fifth iteration step are presented. In Fig. 14 the predicted heat fluxes are displayed. The thermal input of natural-gas was 18.2 MW according to the operating conditions. Due to the pre-heating of the air up to 310 °C an additional input of 1.9 MW was determined. CFD simulation calculated a heat transfer rate to the billets of 10.48 MW which is 52.1% of the total thermal input. This value can be seen as the furnace efficiency. Heat losses of the flue gas and the walls were 5.37 MW and 4.24 MW respectively. But it must be mentioned, that the predicted heat losses at the walls also contains losses due to billet charging and discharging. Fig. 15 shows the temperature and velocity in the furnace along the four burner rows. In the pre-heating zone a gas temperature between 600 and 800 °C was detected and therefore the radiative heat flux to the billets is low in the pre-heating zone with a value between 15 and 39 kW/m2 (see Fig. 16). In the pre-heating zone the convective heat flux reached a value of 3 kW/m2. The maximum convective heat flux in the furnace to the billet occurs between pre-heating and heating zone. A reduction of the cross section due to the water-cooled wall caused an increase of the flue gas velocity which is illustrated in the contour plot of Fig. 15. Higher velocities in this region (between pre-heating and heating zone) constitute a convective heat flux of 5 kW/m2. A comparison with another furnace investigated by Morgado et al. [18] showed that radiative heat fluxes in the pre-heating zone can reach approximately 85 kW/m2 when gas temperature is above 1000 °C in the pre-heating zone. Emadi et al. [5] calculated even higher radiative heat fluxes of about 110 kW/m2 in a walking hearth type furnace with a zone model for the radiative heat transfer. After the billet entered the heating zone the radiative heat flux increased up to 57.5 kW/m2 which is about 95% of the total heat flux in this zone, because the gas temperature was in a range of 1200–1300 °C. Up to the end of the furnace the radiative heat flux reached a maximum of 98% of the total heat flux to the billet. Radiative heat flux in the entire furnace was determined with 93% compared the overall heat transfer at the billet surface. Hence,
Fig. 12. Surface temperature of the last billet for each iteration step.
684
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
Fig. 13. Calculated relative deviations between the iteration steps: (a) Basic vs. first iteration; (b) First vs. second iteration; (c) Second vs. third iteration; (d) Third vs. fourth iteration.
1.6 m). In Fig. 17 the temperature trend calculated by CFD is interrupted at 8.5 m due to the water-cooled wall between pre-heating and heating zone which passes through the furnace at the measurement height of 0.8 m. It was found that a good accordance between measurement and CFD was reached for the entire furnace length. A maximum deviation was detected at the water-cooled wall with 83 K. Basically, the calculated temperatures were too low for all measurement points because of the higher heat losses through the walls including the losses of charging and discharging of the billets. Nevertheless, the mean deviation to the measured temperatures was 42 K. Based on this study it was determined that the CFD model is able to predict the gas phase combustion with high accuracy to the measured values and low computational time due to the use of the steady flamelet approach. Fig. 14. Thermal inputs and heat fluxes in the furnace.
4.4. Heating characteristic of the steel billets it can be concluded that convection had a minor effect on the heating of the billets especially in the heating and soaking zones. In addition to the analysis of the heat fluxes in the furnace, a comparison of the measured and simulated wall temperatures in the furnace was carried out. Measurement was done with the quotient pyrometer for both side walls. The values illustrated in Fig. 17 are the average temperatures of both side walls. The temperature measurement on the wall was done at a height of 0.8 m (overall
Knowledge about the temperature difference inside the billets during the heating is crucial for further processing on the rolling mill. When the billet enters the rolling mill a homogeneous temperature distribution as well as the desired target temperature is important during manufacturing to maintain optimum product quality. Furthermore the maximum temperature difference should be low to avoid deformations of the billets in the furnace. In this
Fig. 15. Contour plots of the furnace: (a) Temperature and (b) Velocity.
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
685
Fig. 16. Convective and radiative heat fluxes on the billets depending on furnace length.
Fig. 17. Measured and calculated wall temperatures in the furnace.
section the temperature distribution in the billet was investigated at different positions in the furnace. The temperature trend was observed along five planes along the furnace. These planes are displayed in Fig. 18 named ‘‘Level A” to ‘‘Level E”. The planes B and D (blue) were placed along the second and third row of burners whereas the other planes (grey) were between the burners. Nine temperature monitors labelled ‘‘t1” to ‘‘t9” were arranged inside the billet for each plane in the CFD model. Also a temperature monitor was set in the middle of the top surface (‘‘t_surf”). The calculated results for ‘‘Level A” only are displayed in Fig. 19, because the nine temperature monitors in all 5 planes showed the same trend along the furnace length. It was found that in the first half of the pre-heating zone the temperature inside the billet is quite homogenous with a small temperature difference. The maximum temperature difference in the billet until the end of the preheating zone increased to 99.9 K between t1 and t6. During the entire process t1 and t6 showed the highest temperature gap. The highest heat fluxes were detected on entering the heating zone and this leads to a faster increase in the temperature, especially for
the top monitors such as t1, t4 and t7. The lowest temperature occurs at position t6 where the heat transfer took place mainly due to conduction from the furnace bottom. The temperature difference reached a global maximum at 12.7 m. At this position the temperature at t6 started to increase rapidly. This is a result of the drop on the specific heat capacity (see Fig. 4) when the temperature is above 830 °C. For this condition the specific heat dropped from over 1000 to 600 J/(kgK) and the position t6 heats up faster than t1. In Fig. 20 the trend of the temperature difference between t1 and t6 is illustrated at all observed planes. Inside the billets the temperature difference started to increase with a constant rate in all planes. At approximately 4 m, when the billets reached the first burners of the pre-heating zone, the difference at the levels A and B were about 15 K lower than the other planes causing a slightly higher heat flux in the middle of the billets. As mentioned above, the temperature difference was determined at approximately 100 K when entering the heating zone. Due to the vast increase in the radiative heat flux, the temperature difference inside the steel increased faster until a maximum of 192 K was reached at
686
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
Fig. 18. Monitor planes through the furnace (grey: between the burners; blue: under the burners). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 19. Temperature trend in the billet at Level A.
Fig. 20. Temperature difference between t1 and t6 for all monitor planes (Level A to E).
level C, whereas a value of 170 K was observed at level B. A first decrease on the temperature difference in the billets was at 13 m which is still part of the heating zone. From this point to the end
of the soaking zone, the temperature differences dropped to a value between 42 and 49 K. Material properties with high dependency on the temperature can be assumed as a crucial part for
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
687
Fig. 21. Billet surface temperature at different positions: (a) 4.9 m; (b) 7.6 m; (c) 10.3 m; (d) 13.8 m; (e) 15.1 m; (f) 17 m.
modelling the heating characteristic and predicting the temperature gradients in the billets. In addition to the temperatures inside the furnace, the temperature trend on the upper surface (labelled with ‘‘t_surf”) was also considered in this study. At the beginning of the pre-heating zone the highest temperature difference occurred at 2.4 m with a value of 25.5 K between the levels A and E due to a higher convective heat flux to the billet. Increasing radiative heat flux in the heating and soaking zone leads to a more homogeneous temperature distribution on the surface. As a consequence the temperature difference becomes lower than 15 K in the heating zone. At the end of the furnace a temperature difference at the billet surface was determined with a maximum value of 7.6 K between level B and C. Contour plots of the surface temperature at different locations inside the furnace are shown in Fig. 21. The position at 4.9 m is below the first burners in the pre-heating zone. At this location CFD predicted higher temperatures at the corners and edges of the billets due to the higher thermal resistivity. A homogeneous surface temperature was determined at the end of the preheating zone due to increasing radiative heat transfer. This is illustrated by Fig. 21(b), which shows the billet at the end of the preheating zone. A higher temperature gradient was detected only at the corners of the billet. The same phenomenon was found at the end of the pre-heating zone at 10.3 m. After the billet was supported to the heating zone the gradient decreases at the edges due to higher exposure on the radiation as it can be seen in Fig. 21(d) at 13.8 m. This billet was located below the first burner row in the heating zone. In the soaking zone a further improvement in the temperature distribution was determined by the CFD simulation. Fig. 21(e) and (f) illustrate the billet at the beginning and the end of the soaking zone where a further decrease on the temperature difference at the surface leads to high product quality. It was pointed out that the higher radiative heat flux had a crucial impact on temperature distribution on the surface, which should be optimised. This could be achieved in future by adapting the furnace for oxygen enriched combustion and an investigation of the heating characteristics of the billets by CFD simulation.
5. Conclusion In this study a novel iterative solution procedure was tested to consider the transient heating process of steel billets in a walking hearth furnace. The iteration was carried out in two steps with a steady-state simulation of the gas phase combustion and the transient heat transfer in the billets. This approach was found to be a good alternative to zone models which calculate the radiative heat flux and convection by assumed zone temperatures and empirical correlations for the heat transfer coefficient (e.g. Kim et al. [4]) or
computationally demanding simultaneous transient simulations of gas phase combustion and heating of the load (e.g. Han et al. [14]). The overall calculation time for the combustion and transient heating was 10 days with consideration of radiative and convective heat transfer. The results of the steady-state simulation predicted that 52.1% (10.48 MW) of the total thermal input was transported to the steel billets. The overall radiative heat transfer to the billets was 93% of these 10.48 MW. The measured wall temperatures inside the furnace showed good accordance to the predicted values with an average deviation of 42 K for the entire furnace length. A constant increase of the temperature difference inside the billets was detected in the pre-heating zone with a value of 99.9 K at the end. The maximum was reached inside the billets in the heating zone where a peak temperature difference of 192 K was reached. This can be explained by the highest heat transfer rate in the furnace which leads to an increased temperature gradient. From this peak value the gradient decreased until the end of the furnace between 42 and 49 K. Optimisation on the furnace operation can be done without expensive test runs. This investigation forms the basis for the future simulation of the furnace under oxygen enriched conditions up to 25 vol% O2 to increase furnace efficiency with low modifications on the furnace. Due to the low enrichment level the product quality is expected to be the same. For this account simulations will be carried out to determine that oxy-fuel combustion can be applied without damage to the furnace or detriment to the desired product quality (temperature distribution). Conflict of interest There is no conflict of interest. Acknowledgment This work was financially supported by the Austrian Research Promotion Agency (FFG), ‘‘Competence-Headquarter: Ökologische und ökonomische Optimierung von industriellen Hochtemperaturprozessen durch Sauerstoffeinsatz” (Project 848281, eCall 4944044). References [1] Intergovernmental Panel on Climate Change (IPPC). IPPC special report on carbon dioxide capture and storage, Cambridge University Press, Cambridge, UK and New York, NY, USA, 2005. [2] A. Jaklic, B. Glogovac, T. Kolenko, B. Zupancic, B. Tezak, A simulation of heat transfer during billet transport, Appl. Therm. Eng. 22 (2002) 873–883. [3] J.H. Jang, D.E. Lee, M.Y. Kim, H.G. Kim, Investigation of the slab heating characteristics in a reheating furnace with the formation and growth of scale on the slab surface, Int. J. Heat Mass Transfer 53 (2010) 4326–4332.
688
R. Prieler et al. / International Journal of Heat and Mass Transfer 92 (2016) 675–688
[4] M.Y. Kim, A heat transfer model for the analysis of transient heating of the slab in a direct-fired walking beam type reheating furnace, Int. J. Heat Mass Transfer 50 (2007) 3740–3748. [5] A. Emadi, A. Saboonchi, M. Taheri, S. Hassanpour, Heating characteristics of billet in a walking hearth type reheating furnace, Appl. Therm. Eng. 63 (2014) 396–405. [6] C.K. Tan, J. Jenkins, J. Ward, J. Broughton, A. Heeley, Zone modelling of the thermal performances of a large-scale bloom reheating furnace, Appl. Therm. Eng. 50 (2013) 1111–1118. [7] A. Jaklic, T. Kolenko, B. Zupancic, The influence of the space between the billets on the productivity of a continuous walking-beam furnace, Appl. Therm. Eng. 25 (2005) 783–795. [8] A.A. Bhuiyan, J. Naser, Numerical modeling of biomass co-combustion with pulverized coal in a small scale furnace, Procedia Eng. 105 (2015) 504–511. [9] A.A. Bhuiyan, J. Naser, Numerical modelling of oxy fuel combustion, the effect of radiative and convective heat transfer and burnout, Fuel 139 (2015) 268–284. [10] A.A. Bhuiyan, J. Naser, Computational modelling of co-firing of biomass with coal under oxy-fuel condition in a small scale furnace, Fuel 143 (2015) 455–466. [11] M.R. Karim, J. Naser, Progress in numerical modelling of packed bed biomass combustion, in: 19th Australasian Fluid Mechanics Conference, Melbourne, Australia, 2014. [12] S.H. Han, S.W. Baek, S.H. Kang, C.Y. Kim, Numerical analysis of heating characteristics of a slab in a bench scale reheating furnace, Int. J. Heat Mass Transfer 50 (2007) 2019–2023. [13] B.F. Magnussen, B.H. Hjertager, On mathematical modeling of turbulent combustion with special emphasis on soot formation and combustion, in: 16th Symposium (Int) on Combustion, The Combust Inst, 1976. [14] S.H. Han, D. Chang, C.Y. Kim, A numerical analysis of slab heating characteristics in a walking beam type reheating furnace, Int. J. Heat Mass Transfer 53 (2010) 3855–3861. [15] S.H. Han, D. Chang, Optimum residence time analysis for a walking beam type reheating furnace, Int. J. Heat Mass Transfer 55 (2012) 4079–4087. [16] M.Y. Gu, G. Chen, X. Liu, C. Wu, H. Chu, Numerical simulation of slab heating process in a regenerative walking beam reheating furnace, Int. J. Heat Mass Transfer 76 (2014) 405–410. [17] V.K. Singh, P. Talukdar, Comparisons of different heat transfer models of a walking beam type reheat furnace, Int. Commun. Heat Mass Transfer 47 (2013) 20–26.
[18] T. Morgado, P.J. Coelho, P. Talukdar, Assessment of uniform temperature assumption in zoning on the numerical simulation of a walking beam reheating furnace, Appl. Therm. Eng. 76 (2015) 496–508. [19] R. Prieler, M. Demuth, D. Spoljaric, C. Hochenauer, Numerical investigation of the steady flamelet approach under different combustion environments, Fuel 140 (2015) 731–743. [20] T.W.J. Peeters, Numerical modeling oft turbulence natural-gas diffusion flames (Ph.D. thesis), Delft Technical University, Delft, The Netherlands, 1995. [21] W. Trinks, Industrial Furnaces, sixth ed., Jon Wiley & Sons, Inc., Hoboken, New Jersey, 2004. [22] A.A. Bhuiyan, A.K.M. Sadrul Islam, M. Ruhul Amin, Numerical study of 3D thermal and hydraulic characteristics of wavy fin and tube heat exchanger, Front. Heat Mass Transfer 3 (2012). [23] A.A. Bhuiyan, M. Ruhul Amin, J. Naser, A.K.M. Sadrul Islam, Effects of geometric parameters for wavy finned-tube heat exchanger in turbulent flow: a CFD modeling, Front. Heat Mass Transfer 6 (2015). [24] A.A. Bhuiyan, M. Ruhul Amin, R. Karim, A.K.M. Sadrul Islam, Plate fin and tube heat exchanger modeling: effects of performance parameters for turbulent flow regime, Int. J. Automot. Mech. Eng. 9 (2014) 1768–1781. [25] T.-H. Shih, W.W. Liou, A. Shabbir, Z. Yang, J. Zhu, A new k-e eddy viscosity model for high reynolds number turbulent flows – model development and validation, Comput. Fluids 24 (3) (1995) 227–238. [26] S.-E. Kim, D. Choudhury, B. Patel, Computations of complex turbulent flows using the commercial code ANSYS FLUENT, in: Proceedings of the ICASE/LaRC/ AFOSR Symposium on Modeling Complex Turbulent Flows, Hampton, Virginia, 1997. [27] ANSYS Fluent Theory Guide 13.0, ANSYS, Inc, Canonsburg, PA, USA, 2010. [28] R. Prieler, M. Demuth, D. Spoljaric, C. Hochenauer, Evaluation of a steady flamelet approach for use in oxy-fuel combustion, Fuel 118 (2014) 55–68. [29] B.E. Launder, G.J. Reece, W. Rodi, Progress in the development of a Reynoldsstress turbulence closure, J. Fluid Mech. 68 (3) (1975) 537–566. [30] N. Peters, Laminar diffusion flamelet models in non-premixed turbulent combustion, Prog. Energy Combust. Sci. 10 (1984) 319–339. [31] N. Peters, Laminar flamelet concepts in turbulent combustion, in: 21st Symposium (Int) on Combustion, The Combust Inst, 1986, pp. 1231–1250. [32] W.P. Jones, J.H. Whitelaw, Calculation methods for reacting turbulent flows: a review, Combust. Flame 48 (1982) 1–26. [33] T.F. Smith, Z.F. Shen, J.N. Friedman, Evaluation of coefficients for the weighted sum of gray gases model, J. Heat Transfer 104 (1982) 602–608.