Analysis of Billet Thermal Behavior and Temperature Setting Optimization in a Walking Beam Reheat Furnace

Abstract

Aiming at the heating process of copper alloy billets, thermal behavior in a walking beam reheat furnace is investigated. A simplified thermal model is presented to predict billet temperature distribution considering specific heat transfer characteristics. Hot gas blackness is corrected basing on composition of chamber combustion atmosphere during the computation of heat flux boundary condition which improves the model accuracy. A fully implicit scheme is employed to solve the discretization of conduction equation in finite volume method. As an inverse application, a practical method to optimize the chamber gas setting temperature is introduced to make a billet temperature track the technical curve well. Design effects are calculated and simulated with Matlab software. Application results for the hot rolled coil production illustrate the effectiveness of the proposed approach.

1. Introduction

The walking beam reheat furnace is a widely used equipment in metallurgical area. Its main role is to heat billets or slabs to a required temperature which is appropriate for hot rolling. It’s a typical longitudinal, multi-phase large furnace which continuously reheats semi-finished billets in metallurgical production line. Technical requirements of the uniform temperature distribution inside a billet are needed at the furnace exit. Generally the maximum temperature difference within a billet should be lower than certain limitation. And a billet temperature need track the setting trajectory well. Those increase the importance of exact prediction for the billet thermal behavior during heating process which determines the product quality greatly.

There already have been some valuable achievements for this type of furnace concerning the heat transfer process. Most of them focus on analyzing thermal characteristics in different heating conditions. The main techniques can be divided into the following two approaches. One is to solve the energy conservation and full Naviere-Stokes equations governing the hot gas flow concerning thermal phenomena in a combustion process.^1,2,3,4,5) Kim et al.¹⁾ introduce a prediction method for transient slab temperature distribution in a walking beam type reheat furnace. Thermal effect is analyzed considering both turbulent reactive flow and radiative heat transfer in the furnace. Hsieh et al.³⁾ perform simulation study by the commercial software STAR-CD aiming at the similar heat transfer problem in a walking beam reheat furnace. Casal et al.⁵⁾ gives a modification design in the energy transport equation and perform a combustion simulation concerning global combustion mechanism that takes the main gaseous species into account. The other technique focuses on the analysis of radiative heat transfer and transient heat conduction equation in a furnace. It usually discusses the thermal behavior by means of model calculation with different dimension governing equations. In some literatures, temperature distribution inside a slab is solved with 1-dimensional thermal model.^6,7,8,9,10) It is a kind of simplified method, but it can reasonably predict thermal behavior suitable to specific combustion process. Galerkin method⁶⁾ and inverse estimation algorithm⁸⁾ for managing the transient term are used. In some other works, the temperature distribution inside a billet is calculated by a 2-dimensional model.^{11,12,13,14,15,16,17)} Basing on this type of model, finite volume method is generally used to solve the radiative heat transfer equation.^11,12) Optimization for the heating schedule and energy conversation are extended and applied.^13,14,15) And comparison research of different thermal models and simplification methods are introduced in some papers.^16,17,18) Meanwhile, a full 3-dimensional temperature field is studied in recent research work.^{19,20,21,22,23,24)} The CFD methodology is combined with 3-dimensional model. Although a CFD method can be used to analyze the furnace thermal characteristics, it really costs much due to large calculation work in solving numerous governing equations and managing complexity of furnace structure as well as the uncertainty of a model. That makes the scheme has certain limits to implement in practical engineering systems considering the accuracy and economy.

From the observation of above literatures, it can be seen that there are various methods available in research of heat transfer process. Each furnace has specific configuration and heating mechanism. There is no uniform scheme suitable to all kinds of reheating furnaces. Thus, considering some constraints in engineering applications, the proper approach needs to be developed aiming at specific thermal characteristics for different furnaces. The present paper focuses on the heating process for copper alloy billets in a walking beam reheat furnace. Prediction of billet thermal behavior and optimal setting method for chamber temperature are presented. Section 2 gives the outline of furnace configuration. The specific burner location results in all heat resource coming from the up side of the furnace. Section 3 introduces the model simplification method and calculation procedure with boundary conditions. Basing on the current thermal model and discretization technique, billet temperature distribution is simulated with Matlab software in Section 4. Meanwhile, a practical method to optimize the chamber temperature setting trajectory is given in Section 5. At last Section 6 gives some conclusions.

2. Configuration of the Walking Beam Reheat Furnace

A walking beam reheat furnace is one of the important heating equipments in metallurgy field. It’s usually used to reheat the cylinder, cuboid alloy billets. The configuration for the furnace considered in the present work is shown in Fig. 1.

Fig. 1.

Configuration of a walking beam reheat furnace (mm).

The furnace is divided into three parts: preheating zone, heating zone and soaking zone. Billets are successively stepped into a furnace where flat-flame burners serve as heating sources. There are totally three-seven burners located on the top of the furnace. There are no burners in preheating zone. Exhaust gas is re-circulated to preheat billets in this part. The heating zone is divided into two sub-zones with eight burners located in each zone. The soaking zone is divided into eight sub-zones with two to four burners in each zone. It’s a logical partition to ensure the temperature uniformity in the chamber during heating process.

There are five fixed skid beams and four walking beams in the present furnace. A fresh billet is transported step by step from inlet to outlet in furnace by walking beams equipped at the furnace bottom. Thus a billet is moved forward by the periodic beam motion. Billet residence time from being charged into the furnace to the exit is about 300 min. After heating process, they reach to the required technical temperature for hot rolling.

3. Mathematical Formulation

3.1. Thermal Characteristics and Conditions

There are intrinsically three types of heat transfer modes including radiation, convection and conduction simultaneously in a heating process. However, because all burners are located on the top of the furnace, no burners being equipped on the furnace walls and bottom, radiation energy emitted by hot gas is absorbed by the upper side of a billet. Then it’s converted into heat energy and further penetrates inside a billet by conduction. And thermal radiation is the dominant heat exchange which counts over 95% heat energy to the billets. Considering above thermal characteristics, the following technical specifications and conditions are given.

(1) Thermal radiation is the only considered mode of heat exchange between the billets and the combustion environment. Only the top surface of the billet serves as interface for heat exchange as all burners are equipped on the furnace top. The heat conduction between billets and beams is negligible, as well as convection heat transfer to the billets.

(2) For computing the radiative heat transfer, the billet width and length are not taken into account. And the temperature of a participating layer is isothermal from the top to bottom along the billet thickness in heating and soaking zones.

(3) As there is no burner in preheating zone, suppose billet temperature is homogeneously distributed when being charged into heating zone after the preheating process.

3.2. Simplification of Heat Conduction Equation

The whole thermal model for a copper billet is composed by two parts: transient heat conduction equation and its boundary conditions. Heat conduction equation mainly describes a billet temperature rises over time. It can be normally calculated with a transient 3-dimensional heat conduction equation as follows:^25,26,27)   

$$\rho c_p\frac{dT}{dt}=\frac{\partial }{\partial x}\left[\lambda \frac{\partial T}{\partial x}\right]+\frac{\partial }{\partial y}\left[\lambda \frac{\partial T}{\partial y}\right]+\frac{\partial }{\partial z}\left[\lambda \frac{\partial T}{\partial z}\right]$$(1)

T represents a billet current temperature and t denotes a specific heat transfer time. x, y, z represent billet width, length and thickness directions respectively. c_p, ρ, λ are billet specific heat, density and conductivity. According to the given conditions in 3.1, burners are all located on the top of the furnace. Therefore the heat energy absorbed by a billet is vertically heat radiation, namely $\frac{\partial T}{\partial x}=0$, $\frac{\partial T}{\partial y}=0$. Therefore, the model can be simplified into 1-dimensional heat conduction equation:   

$$\frac{dT}{dt}=\frac{\lambda }{\rho c_p}\cdot \frac{{\partial }^{2}T}{\partial z^2}=\sigma \mathrm{\hspace{0.17em} }\frac{{\partial }^{2}T}{\partial z^2}\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }$$(2)

where σ is the thermal conductivity coefficient with σ=λ/ρc_p. And (2) describes the billet temperature trend with time along its thickness. It’s the calculating basis for billet temperature prediction in next section 3.4 and 4 after the heat flux boundary condition is computed in section 3.3. Figure 2 illustrates the heat conduction and conversion simplification procedure.

Fig. 2.

Simplification of billet heat conduction.

In Fig. 2, Billet_j represents the jth billet in a furnace, and Δt is a specific period of heating time for a billet to absorb radiation energy. $C_{x,y,z-1}^{i-1},\mathrm{\hspace{0.17em} }\hspace{0.17em}{C}_{x-1,y,z}^{i-1},\mathrm{\hspace{0.17em} }\hspace{0.17em}{C}_{x,y-1,z}^{i-1}$ represent a micro-cell thermal state before heat absorption under 3-dimensional model while $C_{x,y,z}^i$ is the energy state after Δt period of time of heating process. After simplification, the heat transfer only occurs along billet thickness from top to bottom. $C_{z-1}^{i-1}\mathrm{\hspace{0.17em} },C_z^i$ respectively denote the thermal states before and after heat absorption in Δt for a micro-cell with the thickness h under 1-dimensional model.

3.3. Boundary Conditions

The boundary condition needs to manage in order to calculate the temperature distribution inside a billet. Radiative heat transfer is the only considered mode of heat exchange basing on the assumed conditions in 3.1. Especially, high hot gas temperature on the upside of the chamber makes thermal radiation the dominant mode of heat exchange.^6,25) Taking heating zone I for example, all the radiative heat energy to the billet can be given by^25,26)   

$$\mathrm{\hspace{0.17em} }{Q}_b=Q_{t,g\to b}+Q_{t,w\to b}-Q_{bx}$$(3)

where $Q_{t,g\to b}$, $Q_{t,w\to b}$ are separately the billet absorbed radiation energy from hot gas and furnace walls. Q_bx is the billet emitted radiation energy. Thus the final obtained radiation heat flux of the billet is Q_b. By further deducing, we have   

$$Q_{t,g\to b}={\mathrm{E}}_g\cdot S_b$$(4)

  

$$Q_{t,w\to b}=Q_{bx}-Q_{bx}\cdot {\phi }_g$$(5)

where S_b refers to the area sum of all the billet top surfaces in heating zone I as billets are arranged as a row on the beams. φ_g is the angel coefficient between billet top surface and hot gas, also named the shape factor in the chamber. E_g represents the hot gas radiation capability relative to the black body. Then basing on Stefan-Boltzmann law, we have   

$$E_0=C_0\cdot {\left(\frac{T}{100}\right)}^4$$(6)

  

$$E_g={\epsilon }_g\cdot E_0$$(7)

where ε_g is the gas blackness coefficient. E₀ denotes the black body radiation capability with Stefan-Boltzmann constant C₀=5.675[W/(m²·K⁴)]. Then use (5)–(7) to (4) and yield   

$$\mathrm{\hspace{0.17em} }{Q}_b={\epsilon }_g{E}_0{S}_b-Q_{bx}{\phi }_g$$(8)

For the billet emitted radiation energy Q_bx, it can be given by^25,26)   

$$Q_{bx}=\left(\frac{1}{{\phi }_b}-1\right)\cdot Q_b+C_0\cdot \frac{{\epsilon }_b}{{\phi }_b}{\left(\frac{T_b}{100}\right)}^4\cdot S_b$$(9)

where φ_b is the billet angel coefficient to hot gas which indicates the billet radiative absorptivity. T_b is the billet top surface temperature. Using (9) to (8), we have   

$$Q_b=\frac{C_0}{\frac{1}{{\phi }_g}+\frac{1}{{\phi }_b}-1}\cdot \left[\frac{{\epsilon }_g}{{\phi }_g}{\left(\frac{T_g}{100}\right)}^4-\frac{{\epsilon }_b}{{\phi }_b}{\left(\frac{T_b}{100}\right)}^4\right]\cdot S_b$$(10)

By Kirchhoff’s law, we know ε_b and φ_b are consistent under the same temperature.²⁶⁾ If defining the synthetic radiation coefficient ε_s, it can be written as   

$${\epsilon }_s=\frac{1}{\frac{1}{{\phi }_g}+\frac{1}{{\phi }_b}-1}$$

Then we obtain the radiative heat flux equation   

$$Q_b=C_0\cdot {\epsilon }_s\cdot \left[\frac{{\epsilon }_g}{{\phi }_g}{\left(\frac{T_g}{100}\right)}^4-{\left(\frac{T_b}{100}\right)}^4\right]\cdot S_b$$(11)

For ε_g and φ_g, they are not always the same. The angel coefficient φ_g has following relationship with ε_g²⁷⁾   

$$\begin{array}{l}{\phi }_g={\epsilon }_g\left\{1+\frac{S_b/S_f}{1+\left[{\epsilon }_g/(1-{\epsilon }_g)\right]\cdot S_b/S_f}\right\}\\ \mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }={\epsilon }_g\left\{1+\frac{\omega }{1+\left[{\epsilon }_g/(1-{\epsilon }_g)\right]\cdot \omega }\right\}\end{array}$$(12)

S_f is the top and side wall surface area except the furnace bottom area of heating zone I. ω is named as the fender spreading degree with ω=S_b/S_f. For the computation of φ_g, because radiation is not just a local phenomenon, it generally couples the whole dynamical heating environment in a whole chamber. It can take the experienced value used for flame furnaces given in some manuals.^26,27)

Furthermore, ε_g and φ_g also change following by the combustion atmosphere due to the gas radiation characteristics. Thus they need be calculated on basis of the chamber atmosphere. As hot gas radiation has selectivity, the triatomic gases have strong radiation intensity and emissivity. Therefore, the chamber gas composition of CO₂, H₂O plays a key role to the radiation heat flux calculation. The equation for calculating hot gas blackness ε_g is given by^25,27)   

$${\epsilon }_g={\epsilon }_{C{O}_2}+\beta {\epsilon }_{H_{2}O}-\mathrm{\Delta }\epsilon$$(13)

where ${\epsilon }_{C{O}_2},\hspace{0.17em}\hspace{0.17em}{\epsilon }_{H_{2}O}$ represent CO2 and H2O blackness. β denotes the H₂O blackness corrected coefficient. Δε is hot gas blackness corrected coefficient. As known, ${\epsilon }_{H_{2}o},{\epsilon }_{c{o}_2}$ are mainly influenced by hot gas temperature, partial pressure and reservoir thickness. The following procedure is given to calculate the value of hot gas blackness ε_g.

Step1: Calculate the hot gas effective reservoir thickness L.

Basing on the calculation rule for a box space, the chamber volume V_Heat of heating zone I is   

$$V_{Heat}=W_{Inner}\times L_{Heat}\times H_{Heat}$$(14)

where L_Heat, H_Heat and W_Inner represent inner chamber length, height and width of heating zone I. Then the average effective ray distance L can be given by   

$$L=\eta \cdot \sqrt{V_{Heat}}$$(15)

η is the hot gas radiation coefficient with a known value 0.67.²⁵⁾

Step2: Calculate the CO₂, H₂O partial pressure.

According to the combustion situation, suppose CO₂, H₂O partial pressure proportions in the hot gas are known as $P_{C{O}_2},P_{H_{2}O}$. Then the following values can be calculated   

$$P_{H_{2}o}^p=P_{H_{2}O}\cdot L\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} } (m\cdot atm)$$(16)

  

$$P_{C{O}_2}^p=P_{C{O}_2}\cdot L\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} } (m\cdot atm)$$(17)

Using the above calculating results, ${\epsilon }_{C{O}_2},\hspace{0.17em}\hspace{0.17em}{\epsilon }_{H_{2}O}$ can be inquired from the sheet in Industrial Furnace Manual.²⁵⁾

Step3: Calculate corrected coefficient β.

For the corrected coefficient β, it can be inquired from Fig. 3. β changes with $P_{H_{2}O}$ and $P_{H_{2}o}^p$. There are seven curves with different $P_{H_{2}o}^p$ in Fig. 3. We should choose the proper curve basing on the current calculated value of $P_{H_{2}o}^p$. Thus β is determined by chamber combustion state.

Fig. 3.

Corrected coefficient β.

Step4: Calculate hot gas blackness corrected coefficient Δε.

To obtain the corrected coefficient Δε, the following values need be computed first   

$$P_{per}=\frac{P_{H_{2}O}}{P_{H_{2}O}+P_{C{O}_2}}$$(18)

  

$$P_{syn}=(P_{H_{2}O}+P_{C{O}_2})\times L (m\cdot atm)$$(19)

Then Δε can be acquired from Fig. 4.²⁵⁾ Only one group of curves for Δε are given at 520°C. The curves change with temperature rising. Thus Δε is affected by factors including gas temperature, P_syn and P_pers. It should be calculated according to the current thermal state.

Fig. 4.

Corrected coefficient Δε at 520°C.

Step5: With the obtained values, the hot gas blackness ε_g can be calculated by (13). Accordingly φ_g can be obtained by (12).

For the hot gas temperature T_g, shielded thermocouples are equipped on the furnace top to help to calculate its value. Because of the thermal radiation from hot gas to billets and furnace walls, the temperature which is measured by thermocouple is an intermediate value between the local gas temperature and the furnace top temperature in the chamber. To improve the computation accuracy, the following approximate formula is given to calculate the hot gas temperature^6,25)   

$$\mathrm{\hspace{0.17em} }{T}_t^4\approx {\epsilon }_g\mathrm{\hspace{0.17em} }{T}_g^4+(1-{\epsilon }_g\mathrm{\hspace{0.17em} })T_p^4$$(20)

T_t and T_p separately represent the local gas and furnace top temperatures measured by the thermocouples. Accordingly, the hot gas temperature in the furnace can be estimated by   

$$T_g\approx \sqrt[4]{\left[T_t^4-(1-{\epsilon }_g\mathrm{\hspace{0.17em} })T_p^4 \right]/{\epsilon }_g}$$(21)

Using all the above acquired parameters, the billet heat flux density served as boundary condition can be calculated. After the heat flux boundary condition is obtained, the temperature prediction for a billet can be implemented.

3.4. Discretization of Heat Conduction Model

Before predicting temperature distribution inside a billet, the heat transfer equation needs to be discreted both in time and space basing on (2). The discretization is done using a finite volume method. The transient term is managed by the fully implicit scheme. Let the number of divisions of a billet thickness [0, H] be n, where the mesh size along billet thickness is given by h=H/n. Let the number of divisions of the time [0, t_M] be m, t_M represents the heating time. Then the time unit for discretization is given as Δt=t_M/m. Hence we can get the z-t mesh plane shown in Fig. 5.

Fig. 5.

Computational Mesh.

In Fig. 5, T_m,n represents temperature point with coordinate (m·Δt, nh), namely   

$$T_{m,n}=T(m\cdot \mathrm{\Delta }t,nh)$$(22)

Then the temperature points next to T_m,n can be written as   

$$\{\begin{array}{c}{T}_{m-1,n-1}=T[(m-1)\cdot \mathrm{\Delta }t,(n-1)h]\\ T_{m-1,n}=T[(m-1)\cdot \mathrm{\Delta }t,nh] \\ T_{m,n-1}=T[m\cdot \mathrm{\Delta }t,(n-1)h] \\ T_{m+1,n}=T[(m+1)\cdot \mathrm{\Delta }t,nh]\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\\ T_{m,n+1}=T[m\cdot \mathrm{\Delta }t,(n+1)h]\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\\ T_{m+1,n+1}=T[(m+1)\cdot \mathrm{\Delta }t,(n+1)h]\end{array}$$(23)

And the following differentials can be calculated   

$$\{\begin{array}{c}\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\frac{\partial T}{\partial t}|{}_{\begin{array}{l}z=nh\\ t=m\cdot \mathrm{\Delta }t\end{array}}\approx \frac{T_{m+1,n}-T_{m,n}}{\mathrm{\Delta }t},\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\\ \mathrm{\hspace{0.17em} }\frac{\partial T}{\partial z}|{}_{\begin{array}{l}z=nh\\ t=m\cdot \mathrm{\Delta }t\end{array}}\approx \frac{T_{m,n+1}-T_{m,n}}{h},\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\\ \frac{{\partial }^{2}T}{{\partial }^{2}z}|{}_{\begin{array}{l}z=nh\\ t=m\cdot \mathrm{\Delta }t\end{array}}\approx \frac{T_{m,n+1}-2T_{m,n}+T_{m,n-1}}{h^2}\end{array}$$(24)

Use the above (22)–(24) into (2) and yield   

$$T_{m+1,n}=\frac{\sigma \cdot \mathrm{\Delta }t}{h^2}{T}_{m,n+1}+\left(1-2\frac{\sigma \cdot \mathrm{\Delta }t}{h^2}\right)T_{m,n}+\frac{\sigma \cdot \mathrm{\Delta }t}{h^2}{T}_{m,n-1}$$(25)

Then the middle elements in z–t plane can be computed from column 1 to M-1 by difference recurrence formula. The boundary elements are calculated by heat flux and initial boundary conditions. A layer temperature distribution along billet thickness in a period of time t_M can be given by Fig. 6.

Fig. 6.

Layer temperature distribution of a billet.

The computing procedure will be programmed into a function file to operate. In the next section, prediction of a billet temperature distribution basing on above discretization thermal model and boundary condition will be analyzed.

4. Simulation Prediction of a Billet Temperature Distribution

Basing on the above discretization method, simulation analysis of a billet temperature distribution in the heating zone will be discussed with Matlab software. Taking a walking beam reheat furnace project built in Iran for example, technical parameters collected during the running system are given in Table 1.

Table 1. Technological parameters of heating zone I.

Num.	Furnace parameters	Description
1	Burner number	8
2	Burner power	450 kw
3	Heat value of natural gas	38.62 kJ/Nm3
4	Billet material	Iron bronze
5	Beam step velocity	0.065 m/min
6	Billet dimension (Width×Thickness×Length)	330 mm×150 mm×3800 mm
7	Heating zone I dimension (Width×Height×Length)	6200 mm×1600 mm×3800 mm

The billet parameters of ρ, c_p, λ can be inquired from the manuals.^28,29) The current chamber hot gas temperature is measurable. If let the number of divisions along a billet thickness be 100, then the mesh size h=1.5mm. Let the number of division time [0, t_M] be 1000 with the computation time t_M taking as one minute, then the time mesh size Δt=0.06s. Then we can obtain the billet temperature distribution pictures after fifteen minute residence time in heating zone I shown in Fig. 7.

Fig. 7.

Billet temperature prediction after 15 min residence time in heating zone I.

Figure 7(a) illustrates the billet layer temperature distribution along its thickness. Figure 7(b) shows the temperaute rising range at different billet thickness. From the simulation results, it can be found the upper side temperature of a billet is higher than the bottom side. The temperature rising velocity is also faster in the billet top side than the bottom in the current period of time. The average rising velocity along billet thickness is about 13.5°C/min. The temperature difference from the billet bottom to top surface is about 85°C. In order to compare heating effects in different time and position, when the residence time gets to twenty-five minutes, using the difference recurrence formula, simulation results are shown in Fig. 8.

Fig. 8.

Billet temperature prediction after 25 min residence time in heating zone I.

Figure 8 shows temperature gradient reduces to about 45°C in the billet. And the temperature rising velocity at the top side of the billet gets slower than before with the billet staying longer in the furnace. The valid prediction for temperature distribution inside a billet can help to master its temperature rising behavior in different combustion environment so as to improve the billet heat treatment quality in an engineering system.

5. Optimization of Furnace Hot Gas Setting Temperature

Next in turn, a further application of above thermal characteristics will be introduced. Namely calculate the optimal chamber hot gas setting temperature T_g according to the given billet temperature rising curve. For the heat treatment of alloy billets in the reheat furnace, billet temperature is the important parameter to control in each zone, whereas it’s a usual way to accomplish it by controlling the hot gas tempratrue in a practical system. Therefore the gas temperature setting trajectory in the furnace zone has big influence on the quality of heat treatment quality. Basing on the heat transfer characteristics and energy conservation law, we can find the relationship between billet temperature and hot gas temperature. Thereby, the optimization for chamber setting temperature can be achieved, simultaneously reduce the temperature deviation between the required billet temperature and actual values.

5.1. Computation Scheme

From the above analysis of the heat transfer process to a billet, we know the radiation energy is absorbed based on (11) by the billet top surface. Then the radiation energy is changed to heat energy and penetrates from the billet top to bottom by conduction. Gradually the whole billet temperature is raised.

Suppose that the absorbed heat energy is Q_M when billet average temperature is raised from T_i to T_i+1 during a period time of Δτ. Then we have²⁷⁾   

$$\begin{array}{l}{Q}_M=\frac{G}{\mathrm{\Delta }\tau }\cdot \left(c_{p_{i+1}}{T}_{i+1}-c_{p_i}{T}_i\right)\\ \mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }=\frac{G}{\mathrm{\Delta }\tau }\cdot \left(H_{i+1}-H_i\right)\\ \mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }\mathrm{\hspace{0.17em} }=G\cdot \frac{\mathrm{\Delta }H}{\mathrm{\Delta }\tau }\end{array}$$(26)

where G is a billet quality. And H is the billet thermal content at current temperature point. Hence Q_M = Q_b, it can be written as   

$$C_0\cdot {\epsilon }_s\cdot \left[\frac{{\epsilon }_g}{{\phi }_g}{\left(\frac{T_g}{100}\right)}^4-{\left(\frac{T_b}{100}\right)}^4\right]\cdot S_b=G\cdot \frac{\mathrm{\Delta }H}{\mathrm{\Delta }\tau }$$(27)

Basing on (27), we can calculate the temperature setting trajectory for T_g in the case that the billet top surface technical curve is given. The procedure is given as follows.

Step1: Draw up the billet top surface temperature trajectory basing on alloy material and specifications of the current billet.²⁹⁾ This trajectory is also taken as the target temperature curve. For the billet heat treatment time, it can be calculated by the following empirical formula   

$$t_{Heat}=\zeta \cdot H^{1/2}$$(28)

ζ is time factor which has different value with relation to billet material.²⁹⁾ H also represents the billet thickness.

Step2: Calculate the heating delay time of the billet. Then draw up the billet lowest and average temperature curves basing on the delay time. The lowest and average temperature of a billet will keep up with the billet top surface temperature after the delay time τ_L and τ_A respectively. The time delay changes with different billet sizes and heating modes. The calculating formula of delay times for several widely used types of flame furnaces can be acquired from industrial manuals.²⁷⁾

Step3: Divide the total heating time into N equal amounts. Accordingly the time points of t₁, t₂, …, t_N are obtained.

Step4: Acquire the average temperature T_a(1, 2) between t₁ and t₂ from the billet average temperature trajectory. Calculate the billet average thermal content with G, ΔH, Δτ at the current temperature basing on (26).

Step5: Then basing on (27), the chamber average temperature T_f (1, 2) in Fig. 9, namely the hot gas temperature T_g can be calculated.

Fig. 9.

Illustration of furnace hot gas setting temperature.

Step6: Repeat above procedures from i= 1 to i= N. And furnace hot gas temperature setting trajectory can be drawn up through connecting the points.

At last some corrections for the setting trajectory are needed according to expert experience, system simulation and running effects. The procedure can be illustrated with Fig. 9.

For the four trajectories in Fig. 9, the target temperature curve also known as the billet top surface temperature trajectory is given first. Then the lowest temperature curve can be obtained. Here the lowest temperature is taken as the billet bottom temperature. Next, the average temperature trajectory will be calculated. At last the hot gas temperature setting trajectory in the furnace can be drawn up. The main purpose for the optimization of hot gas setting temperature is to minimize the temperature difference inside a billet during heating process .

5.2. Temperature Setting Optimization in Heating Zone I

The hot gas temperature setting trajectory in a furnace has big influence on the billet temperature rising trend and product quality. Take the heating process in heating zone I into account. The computation is done considering a billet made of beryllium bronze with the size 380 mm×205 mm×2510 mm(Width×Thickness×Length) being processed. Use the above computing scheme. Figure 10 shows the corresponding temperature trajectories after calculation with Matlab software.

Fig. 10.

Furnace hot gas temperature setting optimization.

There are seven trajectories in Fig. 10. Sampling period is taken as 1.5 min. For the conventional setting temperature method, it’s a fixed value shown as Former set in Fig. 10. While using the new computation method, the hot gas setting temperatures is calculated in program according to the Target top curve under the current heating situation. The Target top is known as a billet top surface temperature trajectory which is given first according to the billet material and specificatoin as stated in step1.²⁹⁾ Basing on the Target top curve, New bottom and Average trajectories can be obtained according to the computing procedure step2. Then New set, namely the hot gas temerature setting trajectory in the furnace can be determined. To compare with the former setting method, we deduce the billet top surface and bottom temperature trajectories in reverse supposing that Former set is taken as the hot gas temperature setting curve. And then Former top and Former bottom trajectories are calculated as shown in Fig. 10. From the picture, it can be noticed that the values in Former top are higher than those in Target top curve. Moreover, the temperature deviations between Former top and Former bottom are larger than those calculated by Target top and New bottom along the billet thickness. Especially, during the first period of time in the heating process, the temperature deviation obtained by former method is more obviously big than that obtained by the new scheme. This can be demonstrated futher by the effects of billet heat treatment in the next part 5.3. It’s a practical method applicable to amend the furnace gas temperature setting trajectory when the billet type or heating environment is changed.

5.3. Application Effects

To illustrate the application effects of the proposed method, the comparison parameters of heat treatment results and product mechanical properties are presented. As known, alloy billets are heated to the processing temperature through sequential furnace zones. Then they are fed to the mill for hot rolling. As one of important technical parameters, the billet section temperature difference is required less than 5°C at the furnace exit. And the chamber gas setting temperature curve plays a key role on it. Take several batches of beryllium bronze blocks for heat treatment.The alloy billets have the same thickness of 205 mm. Basing on the statistical data from the site engineers recorded in the project, the processing parameters are listed in Tables 2 and 3.

Table 2. Billet parameters by former temperature setting method.

Billet size (Width×Thickness×Length)	Weight (Kg)	Tapping temperature	Section temperature difference
380 mm×205 mm×2515 mm	2465	802°C	4.5°C
375 mm×205 mm×2485 mm	2405	782°C	4.8°C

Table 3. Billet parameters by new temperature setting method.

Billet size (Width×Thickness×Length)	Weight (Kg)	Tapping temperature	Section temperature difference
425 mm×205 mm×2690 mm	2535	792°C	3.6°C
380 mm×205 mm×4365 mm	2690	789°C	3.3°C

Table 2 gives the heat treatment parameters for two batches of alloy billets when produced by the former chamber gas temperature setting method. The temperature setting value is 850°C in heating zone and 820°C in soaking zone. The statistical billet section temperature differences at the furnace exit are separately 4.5°C and 4.8°C. Table 3 presents billet heating parameters when using the new temperature setting method. The average temperature deviations in billets are separately 3.6°C and 3.3°C at furnace exit. Although all the final parameters can meet technical requirements, namely billet temperature gradient no more than 5°C, deviation values in Table 3 are obviously decreased compared with those in Table 2.

Moreover, because billet thermal stress and crystalline state are also affected by its temperature rising process during heat treatment, the chamber hot gas temperature setting curve also acts on the mechanical properties of hot rolled coil products. The testing mechanical properties for the coils produced using the former temperature setting method in the production line are given in Table 4. In view of the possible influence that different test-pieces bring to the results, it implements the standard ISO6892.1-2009 for the mechanical property test.³⁰⁾

Table 4. Coil mechanical properties by former temperature setting method.

Nominal thickness (mm)	Tensile strength 0b (MPa)	Yield strength 0.2% (MPa)	Elongation (%)
12	≥ 680	≥ 645	≥ 5.5

When the new chamber temperature setting method is applied, corresponding mechanical properties of strip coil products with the same nominal thickness are given in Table 5.

Table 5. Coil mechanical properties by new temperature setting method.

Nominal thickness (mm)	Tensile strength 0b (MPa)	Yield strength 0.2% (MPa)	Elongation (%)
12	≥ 715	≥ 660	≥ 5.7

From the parameters in Tables 4 and 5, it can be seen that mechanical properties are improved when using the new temperature setting trajectory. The strip coil tensile strength is enhanced about 5.1% than the former method. Yield strength and elongation also have been increased at least 2.3% than before. Furthermore the fuel consumption can be saved when using the new method in applications. And the exact hot gas temperature curve also helps to maintain the billet temperature rises smoothly which can achieve good surface glossiness for strip coil products.

6. Conclusion

A simplified methodology is used to investigate the thermal characteristics in a walking beam reheat furnace concerning the specific heating mechanism. A practical method is introduced to set the furnace hot gas temperature. The Matlab software is used to simulate and compute the design effects. Actural application performances are given and analyzed. In conclusion, the following results are summarized:

(1) Burner location determines the thermal exchange mode and heat flow direction in the present reheat furnace. It also determines the calculating scheme of heat transfer model.

(2) The thermal radiation is the main heat exchange mode in the current reheat furnace. The chamber gas blackness ε_g has influence on the radiation heat flux calculation. The gas composition should be analyzed first as triatomic gases have strong radiation intensity in high temperature.

(3) Billet temperature distribution is calculated by using a finite volume discretization method. The simulation effect shows temperature rising velocity on the up side is faster than the bottom part inside a billet during the first period of time. With the decreasing of the billet surface temperature rising velocity, the billet central layer temperature continues to rise and the billet temperature gradient reduces gradually.

(4) The chamber hot gas setting temperatures can be computed when expected billet top surface temperature rising trajectory and billet specifications are given. The method is practical which can implement conveniently in an engineering system.

Acknowledgments

This work was supported by the National Science of China (grant no. 61401055), funding of Liaoning Science Technology Research Program of China (grant no. L2013463). The authors will also express sincere thanks to professor Gwynfor Richards of Brandon University and senior engineer Xing Haibin who give much good advice during the preparation of the manuscript which improves the presentation.

References

• 1)   J. G.  Kim and  K. Y.  Huh: ISIJ Int., 40 (2000), 1115.
• 2)   C. T.  Hsieh,  M. J.  Huang, S. T. Lee and C. H. Wang: Numer. Heat Transf. A, 57 (2010), 1.
• 3)   C. T.  Hsieh,  M. J.  Huang and  S. T.  Lee: Numer. Heat Transf. A, 53 (2008), 966.
• 4)   J. G.  Kim,  K. Y.  Huh and  I. T.  Kim: Numer. Heat Transf. A, 38 (2000), 589.
• 5)   J. M.  Casal,  J.  Porteiro,  J. L.  Míguez and  A.  Vázquez: Appl. Therm. Eng., 86 (2015), 69.
• 6)   A.  Steinboeck,  D.  Wild and  T.  Kiefer: Int. J. Heat Mass Transf., 53 (2010), 5933.
• 7)   J.  Harish and  P.  Dutta: Ironmaking Steelmaking, 32 (2005), 151.
• 8)   C.  Miao,  Y.  Kai and  Y. F.  Liu: J. Iron Steel Res. Int., 19 (2012), 13.
• 9)   M. Y.  Gu,  G.  Chen and  X. H.  Liu: Int. J. Heat Mass Transf., 76 (2014), 405.
• 10)   A.  Steinboeck,  D.  Wild,  T.  Kiefer and  A.  Kugi: Math. Comput. Simul., 82 (2011), 392.
• 11)   M. Y.  Kim: Int. J. Heat Mass Transf., 50 (2007), 3740.
• 12)   S. H.  Han,  D. J.  Chang and  C.  Huh: Energy, 36 (2011), 1265.
• 13)   M.  Suzuki,  K.  Katsuki and  J. I.  Imura: J. Process Control., 24 (2014), 225.
• 14)   A.  Saboonchi,  S.  Hassanpour and  S.  Abbasi: Energy Convers. Manag., 49 (2008), 3211.
• 15)   J. Y.  Jang and  J. B.  Huang: Appl. Therm. Eng., 85 (2015), 313.
• 16)   V. K.  Singh and  P.  Talukdar: Int. Commun. Heat Mass Transf., 47 (2013), 20.
• 17)   N.  Depree,  J.  Sneyd and  S.  Taylor: Comput. Chem. Eng., 34 (2010), 1849.
• 18)   J. H.  Jang,  D. E.  Lee and  M. Y.  Kim: Int. J. Heat Mass Transf., 53 (2010), 4326.
• 19)   A.  Emadi,  A.  Saboonchi and  M.  Taheri: Appl. Therm. Eng., 63 (2014), 396.
• 20)   M. H.  Bordbar,  K.  Myohanen and  T.  Hyppanen: Appl. Therm. Eng., 76 (2015), 344.
• 21)   J. M.  Casal,  J.  Porteiro and  J. L.  Miguez: Appl. Therm. Eng., 86 (2015), 69.
• 22)   S. H.  Han,  S. W.  Baek and  M. Y.  Kim: Int. J. Heat Mass Transf., 52 (2009), 1005.
• 23)   D. N.  Trivic and  C. H.  Amon: Int. J. Heat Mass Transf., 51 (2008), 2711.
• 24)   S. H.  Han and  D. Y.  Chang: Int. J. Heat Mass Transf., 55 (2012), 4079.
• 25)   B. Q.  Wang: Industrial Furnace Manual, 3rd Ed., China Machine Press, Beijing, (2012), 23.
• 26)   Q. F.  Cai: Heating Furnace, Metallurgical Industry Press, Beijing, (1983), 35.
• 27)   B. Q.  Wang: Industrial Furnace and Energy-saving Manual, Chemical Industry Press, Beijing, (2006), 62.
• 28)   X. C.  Xu and  L. X.  Zhou: Combustion Technology Manual, Chemical Industry Press, Beijing, (2007), 624.
• 29)   W. J.  Zhong: Copper Fabrication Manual, Metallurgical Industry Press, Beijing, (2009), 125.
• 30)  ISO6892.1: 2009, Metallic Materials-Tensile Testing-part 1: Method of Test at Room Temperature.