Temperature and Thermal Stress Analysis of a Hot Blast Stove with an Internal Combustion Chamber

Abstract

In this study, the temperature and thermal stress fields of an internal combustion hot blast stove were calculated and analysed. Turbulent, species transport, chemical reaction, radiation, and porous media models were implemented in a computational fluid dynamics model. Thermal boundary conditions on the structure of the hot blast stove were calculated based on the analytic adiabatic Y-plus method. A method to interpolate the thermal boundary conditions to a finite element mesh was developed, and the boundary conditions were mapped through the proposed method. In the on-gas period, the vortex was generated in the dome, and it made the variation of the temperature field in the checker chamber. The maximum temperature of the flue gas reached 1841 K in the on-gas period. In the on-blast period, the flow was considerably even compared to the on-gas period, and the average blast temperature reached 1345 K. The outer region of the checker chamber is shown to be continuously exposed to a higher temperature, which makes the region the main domain in managing the deterioration of the refractory linings. The shell temperature did not change during the operation due to the lower thermal diffusivity of the refractory linings, where the inner surface of the refractory had a maximum temperature change from 1441 K to 1659 K. The maximum temperature of the shell was 418.4 K at the conical region of the checker chamber side. The conical region had the higher maximum and middle principal thermal stresses due to the presence of a large temperature gradient around the conical region, where the largest maximum and middle principal stresses were 300.6 MPa and 192.0 MPa, respectively. The conical region was found to be a significant area of interest where it had a higher temperature and thermal stress.

1. Introduction

A hot blast stove, a facility that blows preheated air into a blast furnace, has a vital role in the ironmaking process. The types of hot blast stoves are separated as internal, external, and top combustion types in accordance with a location and shape of a combustion chamber. The internal combustion type has combustion and checker chambers separated by a partition wall in a single tower. The external type has combustion and checker chambers in each separated tower connected from their top. The top combustion type has the combustion and checker chambers in a single tower, but the combustion chamber locates in the dome. The operating cycle of the hot blast stove consists of on-gas and on-blast periods. In the on-gas period, fuel and air come into a burner and combust, and the generated flue gas is blown to the checker chamber to store the heat in the checker bricks. In the on-blast period, a cold blast is blown to the checker chamber, absorbs heat from the checker bricks, and passes to the blast furnace through the hot blast outlet. Due to this long-term periodic operation, the thermal gradient and stress are generated on the hot blast stove. In addition, the thermal stress causes stress corrosion cracking along with the corrosive atmosphere created by the flue gas. The stress corrosion cracking can be evaluated using non-destructive testing methods such as radiographic or ultrasonic tests [1]. However, it is impractical to inspect all regions since the scale of the hot blast stove is enormous. In this circumstance, it is crucial to limit the area to inspect by analysing which part would have the higher temperature and the higher tensile principal stress. In addition, the thermal stress analysis is important because there are reports of explosions due to the stress corrosion cracking in the hot blast stove applications [2,3].

The performance, fluid flow, and temperature of the hot blast stove have been studied in numerous studies. Zhang et al. [4] analysed the combustion, gas flow, and heat transfer processes for the huge-sized dome combustion hot blast stove, which was designed to have a blast temperature of 1300 °C. The process flow and layout were optimised through the numerical simulations, and the average blast temperature reached 1300 °C with single blast furnace gas (BFG) thanks to the developed numerical model. Rieger et al. [5] successfully modelled the fuelling and combustion process of blast furnace gas with the eddy dissipation concept model considering the reaction mechanism of CH₄/air consisting of 41 elementary reactions. They found that the peak temperature within the dome was an important process control parameter to reduce pollutant formation. Guo et al. [6] studied the effect of the number of burner nozzles on the flow and temperature fields in the combustion chamber. It is shown that an adequate number of nozzles in the combustion chamber will give the best performance from their study. Qi et al. [7] utilised a 3D numerical model to optimise the structure of the annual cavity ring of the top combustion hot blast stove. The proposed structural design reduced an asymmetric swirling flow. It increased the resultant highest temperature by around 40 K. Zetterholm et al. [8] proposed a dynamic model for the hot blast stove operation and the concept of stove oxygen enrichment with flue gas recirculation (SOE-FGR) that operates the blast stove without any coke oven gas (COG). The dynamic model analysed the proposed concept, and the concept was shown to increase the blast temperature while improving the efficiency of the entire process. Zhang et al. [9] developed a 3D numerical model to simulate the fluid flow of the top combustion hot blast stove. The design of the neck was modified to resolve the attachment phenomenon of the fluid that occurred in the stove, and the effectiveness of the modified design was verified via numerical simulations.

The stress analysis is also one of the issues in the hot blast stove studies. Yan and Cheng [10] analysed the stress of the top combustion hot blast stove shell. Their numerical model assumed that the fluid’s pressure and the refractory lining’s thermal expansion are constant, making the analysis imprecise. Hwang and Lee [11] studied about temperature and thermal stress of the refractories of the external combustion hot blast stove. They computed thermal stress considering the mortar filled between the refractories and studied the change in the thermal stress according to the number of linings and the existence of the mortar. Gan et al. [12] analysed the temperature and thermal stress of the refractories and the shell of the external combustion hot blast stove. A 3D model was developed to calculate the temperature and velocity fields during the operation. Their model assumes that the surface temperature of the hot blast stove structure and the adjacent fluid temperature are identical.

The hot blast stove can be categorised as an industrial pressure vessel [13]. Understanding the impact of thermal stresses on pressure vessels is crucial, as it is closely linked to hot blast stoves. Several studies have used finite element analysis (FEA) to investigate thermal stresses in pressure vessels, focusing on various factors such as vessel material, geometry, and operating conditions. Rajić et al. [14] conducted a study to optimise a hot water fire-tube boiler using finite element analyses. The kriging response surface was constructed to determine the correlation between the stresses of the tube plate and its thickness. The thickness of the tube was found to be the most critical for the optimisation. Jaremkiewicz et al. [15] proposed a new design thermometer to measure transient fluid temperatures in thick-walled pressure vessels and a method to predict thermal stresses with a finite element model. The proposed design was able to measure the transient fluid temperature with adequate accuracy and be adapted to a wide range of objects. They showed that temperature measurement accuracy was crucial for stress calculations, especially for rapid temperature changes. Taler et al. [16] presented two methods for monitoring thermal stresses in thick-walled pressure components of steam boilers. The first method involves the inverse heat conduction problem to determine the inner surface temperature from the measured outer insulated surface temperature. It does not need the temperature of the fluid and the heat transfer coefficient. The second method uses a numerical-experimental method to determine the transient heat transfer coefficient with the 3D-inverse heat conduction problem. The methods were applied to the thick-walled horizontal steam header to validate the proposed methods. Hoseinzadeh and Heyns [17] analysed the thermo-structural fatigue lifetime of a heat exchange based on CFD and FEA with ASME VIII Boiler and Pressure Vessel Code. They showed that the tubes connected to the tube sheet were the critical points that have the shortest lifetime. The equivalent thermal stress was determined, and its corresponding lifetime was evaluated for each loading condition. Rajić et al. [18] presented experimental and numerical analyses of stress and strain of boiler tube plate. The numerical results were verified with the strain measurement by 3D digital image correlation, not only with the conventional strain measurement using strain gauges. Yıldırım et al. [19] presented a nonlinear model to calculate the thermal stress of functionally graded spherical pressure vessels. The Halpin–Tsai homogenisation scheme was adopted in their study to describe temperature-dependent material properties of the graded material. They combined Pseduspectral Chebyshev Method and the fixed-point method to analyse the thermal stress, and the result was compared and verified with a commercial finite element program. Jeong et al. [20] proposed simple prediction equations of the stress intensity factor via the thermal load under cool-down and heat-up conditions for reactor pressure vessels. They conducted 360 cases of finite element analyses for the different nozzle shapes and heating/cooling conditions. The stress intensity factor solutions were obtained using curve fitting of the finite element analyses solutions. The proposed equations showed good agreement with the finite element analysis solutions.

The thermal stress analysis has yet to be conducted for internal combustion hot blast stoves, while the thermal stresses of the dome and external combustion hot blast stoves have been studied. In addition, the studies on the fluid flow and temperature analyses mainly focused on the dome and external combustion hot blast stoves. In this study, to investigate the areas of interest for managing refractory linings and the shell, a numerical model was established to analyse the flow phenomenon, structure temperature, and thermal stress during an internal combustion hot blast stove operation. The fluid flow was analysed by means of the computational fluid dynamics (CFD) simulation. Turbulent, species transport, chemical reaction, radiation, and porous media models were implemented in the CFD model. Thermal boundary conditions on the structure of the hot blast stove were calculated with the CFD result based on the analytic adiabatic Y-plus method [21]. Calculated thermal boundary conditions were mapped through the developed mapping method, which can interpolate the value at a cell centre of the CFD grid to a finite element analysis (FEA) mesh node. Structure temperature and thermal stress were calculated by the finite element method (FEM). The flow characteristics and regions with high thermal stresses were investigated based on the numerical result. The significant areas of interest for refractory linings and the shell were pointed out from the investigations of the numerical result.

2. Numerical Modelling

The internal combustion hot blast stove has combustion and checker chambers and a dome, as shown in Figure 1a. The conical region is a part that has a conical shape that connects the dome and the tower. The hot blast stove has refractory linings, and it is encased in a metal shell. In the on-gas period, the air and fuel gases come through the air and gas inlets which are in the bottom of the combustion chamber, and the combusted hot flue gas store its heat into the checker bricks and goes out through the flue gas outlet, as shown in Figure 1b. In the on-blast period, cold air comes through the cold air inlet, which is the same as the flue gas outlet in the on-gas period, and the hot blast is generated by the cold air absorbing the heat from the checker bricks and goes out through the hot blast outlet, as shown in Figure 1c. The hot blast stove is heated due to the flue gas and the hot blast during the operation, which generates thermal stress on the shell.

A numerical model was developed to analyse the flow phenomenon, temperature, and thermal stress fields of the hot blast stove during the operation. The 3D fluid flow and temperature in the hot blast stove were analysed with the commercial software Ansys Fluent. The thermal boundary conditions on the structure are calculated based on the CFD results and mapped to the FEA mesh. The temperature and thermal stress fields of the structure were analysed by the finite element method with the commercial software Abaqus.

2.1. CFD Analysis of a Hot Blast Stove

2.1.1. Fluid Flow Model

In the fluid flow in the hot blast stove, the following assumptions are established: (a) the fluid is Newtonian fluid, and (b) it follows an incompressible ideal gas law, and (c) turbulent flow regime. The incompressible ideal gas model is known to be valid with a Mach number less than 0.3 [22]. Since the flow in the hot blast stove adopted in this study did not have a Mach number greater than 0.3, the incompressible ideal gas model was implemented, while a fully compressible model requires high computational cost. In addition, it has been implemented for various studies on the flow that involves a combustion process [23,24,25]. With the assumptions, the continuity, linear momentum, and energy equations are described as follows:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \rho u_i}{\partial x_i}=0$$

(1)

$$\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right]+\rho g_i$$

(2)

$$\begin{array}{c}\frac{\partial \left(\rho c_{p}T\right)}{\partial t}+\frac{\partial \left(\rho c_p{u}_{i}T\right)}{\partial x_i}=\frac{\partial }{\partial x_i}\left(\lambda \frac{\partial T}{\partial x_i}-\rho c_p\overline{u_i^{\prime }{T}^{\prime }}\right)\end{array}$$

(3)

where t is time, $\rho$ is the fluid density, $u_i$ is the velocity component, $x_i$ is the Cartesian coordinate component, p is the pressure, $\mu$ is the dynamic viscosity, $g_i$ is the gravity component, T is the fluid temperature, $c_p$ is the fluid specific capacity, and $\lambda$ is the fluid thermal conductivity. The standard k-$\epsilon$ model was adopted to simulate the turbulent flow of the hot blast stove. Since the hot blast stove has no recirculation or swirling flows, the standard k-$\epsilon$ model was regarded as appropriate for turbulence modelling. It has been implemented in numerous similar applications [5,12,26,27]. The transport equations to compute the turbulent kinetic energy k and the turbulent kinetic energy dissipation $\epsilon$ are described as follows [28]:

$$\begin{array}{c}\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho k{u}_i\right)}{\partial x_i}=\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]+G_k+G_b-\rho \epsilon \end{array}$$

(4)

$$\begin{array}{c}\frac{\partial \left(\rho \epsilon \right)}{\partial t}+\frac{\partial \left(\rho \epsilon u_i\right)}{\partial x_i}=\frac{\partial }{\partial x_i}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_i}\right]+C_{1\epsilon }\frac{\epsilon }{k}(G_k+C_{3\epsilon }{G}_b)-C_{2\epsilon }\rho \frac{{\epsilon }^2}{k}\end{array}$$

(5)

where $C_{1\epsilon }$, $C_{2\epsilon }$, and $C_{3\epsilon }$ are the model constants. ${\sigma }_k$ and ${\sigma }_{\epsilon }$ are the turbulent Prandtl numbers for k and $\epsilon$. ${\mu }_t$ is the turbulent viscosity, respectively. $G_k$ and $G_b$ are the generation of turbulent kinetic energy due to the mean velocity and buoyancy, respectively. $C_{3\epsilon }$ is the constant that determines the effect of the buoyancy to $\epsilon$. ${\mu }_t$, $G_k$, $G_b$ and $C_{3\epsilon }$ are defined as follows:

$$\begin{array}{cc}\hfill {\mu }_t& =\rho C_{\mu }\frac{k^2}{\epsilon }\hfill \\ \hfill G_k& =-\rho \overline{u_i^{\prime }{u}_j^{\prime }}\frac{\partial u_j}{\partial x_i}\hfill \\ \hfill G_b& =-g_i\frac{{\mu }_t}{\rho {\mathrm{Pr}}_t\frac{\partial \rho }{\partial x_i}}\hfill \\ \hfill C_{3\epsilon }& =tanh\left|u_g/u_{gp}\right|\hfill \end{array}$$

(6)

where $C_{\mu }$ is the model constant, and ${\mathrm{Pr}}_t$ is the turbulent Prandtl number for energy. $u_g$ and $u_{gp}$ are the velocity components parallel and perpendicular to the gravitational vector, respectively. The k-$\epsilon$ model parameters were specified as follows: $C_{1\epsilon }=1.44$, $C_{2\epsilon }=1.92$, $C_{\mu }=0.09$, ${\sigma }_k=1.0$, ${\mathrm{Pr}}_t=0.85$.

2.1.2. Species Transport and Chemical Reaction Models

The species transport model was adopted to analyse the mixing characteristics of the fuel gas and air. The local mass fraction of the species i, $Y_i$, is computed by solving the following convection–diffusion equation [29]:

$$\frac{\partial \left(\rho Y_i\right)}{\partial t}+\nabla ·\left(\rho \mathbf{u}{Y}_i\right)=-\nabla ·{\mathbf{J}}_i+R_i$$

(7)

where u is the velocity vector, ${\mathbf{J}}_i$ is the diffusion flux, and $R_i$ is the net rate of production of species i. The diffusion flux is derived due to the gradients of concentration and temperature, and it can be written as Equation (8) using Fick’s law.

$${\mathbf{J}}_i=-\left(\rho D_{i,m}+\frac{{\mu }_t}{{\mathrm{Sc}}_t}\right)\nabla Y_i-D_{T,i}\frac{\nabla T}{T}$$

(8)

$D_{i,m}$ is the mass diffusion coefficient for species i in the mixture, $D_{T,i}$ is the thermal diffusion coefficient, and ${\mathrm{Sc}}_t$ is the turbulent Schmidt number specified as 0.7 in this study.

The eddy-dissipation model [30] was employed to handle the combustion process. The eddy-dissipation model assumes that the combustion is mixing-limited, which means the rate of the chemical reactions is relatively more rapid than the rate of turbulent mixing. $R_{i,r}$, the net rate of production of species i due to reaction r, is given by the smaller value of the values computed with Equations (9) and (10).

$$\begin{array}{c}\hfill R_{i,r}=v_{i,r}^{\prime }{M}_{w,i}A\rho \frac{\epsilon }{k}\underset{R}{min}\frac{Y_R}{v_{R,r}^{\prime }{M}_{w,R}}\end{array}$$

(9)

$$\begin{array}{c}\hfill R_{i,r}=v_{i,r}^{\prime }{M}_{w,i}AB\rho \frac{\epsilon }{k}\frac{{\sum }_P{Y}_P}{{\sum }_j^N{v}_{j,r}^{″}{M}_{w,j}}\end{array}$$

(10)

$Y_P$ is the mass fraction of product species P, $Y_R$ is the mass fraction of a particular reactant R, $M_{w,i}$ is the molecular weight of species i, $v_{i,r}^{\prime }$ is the stoichiometric coefficient for reactant i, and $v_{i,r}^{″}$ is the stoichiometric coefficient for a product i in reaction r. A and B are empirical constants equal to 4.0 and 5.0, respectively. The large eddy mixing time scale $k/\epsilon$ is a dominant factor of the chemical reaction in the model. The combustion occurs when turbulence is present ($k/\epsilon >0$), and an ignition source is not required to initiate it.

2.1.3. Radiation Model

Radiation heat transfer has a significant role in heat exchange processes. The radiation was implemented with the discrete ordinates (DO) radiation model. The DO radiation model has been widely used for numerous engineering applications [31]. It provides a good compromise between accuracy and computational cost [32] and is applicable for complicated geometries [33,34]. The radiative transfer Equation (RTE) of the DO model is expressed as follows [35]:

$$\nabla ·\left(I(\mathbf{r},\mathbf{s})\mathbf{s}\right)+(a+{\sigma }_s)I(\mathbf{r},\mathbf{s})=a{n}^2\frac{\sigma T^4}{\pi }+\frac{{\sigma }_s}{4\pi }{\int }_0^{4\pi }I(\mathbf{r},{\mathbf{s}}^{\prime })\mathsf{\Phi }\left(\mathbf{s}·{\mathbf{s}}^{\prime }\right)\mathbf{d}{\mathsf{\Omega }}^{\prime }$$

(11)

where I is the radiation intensity that depends on the position r and direction s, ${\mathbf{s}}^{\prime }$ is the scattering direction vector, a is the absorption coefficient, n is the refractive index, ${\sigma }_s$ is the scattering coefficient, $\mathsf{\Phi }$ is the phase function, ${\mathsf{\Omega }}^{\prime }$ is the solid angle, and $\sigma$ is Stefan-Boltzmann constant.

2.1.4. Porous Media Model

The checker bricks in the checker chamber were considered non-equilibrium porous media to simplify the model, and several studies were conducted based on this assumption [7,9,12]. For single-phase flow through a porous media, the volume-averaged mass and momentum conservation equations are described as follows [36]:

$$\begin{array}{c}\frac{\partial \left(\gamma \rho \right)}{\partial t}+\nabla ·\left(\gamma \rho \mathbf{u}\right)=0\end{array}$$

(12)

$$\begin{array}{c}\begin{array}{cc}\hfill \frac{\partial \left(\gamma \rho u_i\right)}{\partial t}+\frac{\partial (\gamma \rho u_i{u}_j)}{\partial x_j}=& -\frac{\partial \left(\gamma p\right)}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\gamma \mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\gamma \rho \overline{u_i^{\prime }{u}_j^{\prime }}\right]\hfill \\ & +\gamma \rho g_i-\left(\frac{{\gamma }^2\mu }{K}+\frac{{\gamma }^3{C}_2}{2}\rho \left|\mathbf{u}\right|\right)u_i\hfill \end{array}\end{array}$$

(13)

where $\gamma$ is the porosity, K is the permeability, and $C_2$ is the inertial resistance factor. The energy equations of the fluid and solids are described as follows:

$$\begin{array}{cc}\hfill \frac{\partial \left(\gamma \rho c_{p}T\right)}{\partial t}+\frac{\partial \gamma \rho c_p{u}_{i}T}{\partial x_i}& =\frac{\partial }{\partial x_i}\left(\gamma \lambda \frac{\partial T}{\partial x_i}-\gamma \rho c_p\overline{u_i^{\prime }{T}^{\prime }}\right)+h\beta (T_s-T)\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \frac{\partial \left[(1-\gamma ){\rho }_s{c}_{p,s}{T}_s\right]}{\partial t}& =\frac{\partial }{\partial x_i}(1-\gamma ){\lambda }_s\frac{\partial T_s}{\partial x_i}+h\beta (T-T_s)\hfill \end{array}$$

(15)

where subscript s denotes the solid, $\beta$ is the surface area density, and h is the heat transfer coefficient for the interface between the fluid and the solid. The shape of the checker brick used in this study is represented in Figure 2. The hole diameter D is 52 mm, and the distance between the holes L is 78 mm. The porosity and the surface area density were calculated with the following equations.

$$\begin{array}{cc}\hfill \gamma & =\frac{\pi }{2\sqrt{3}}\frac{D^2}{L^2}\hfill \\ \hfill \beta & =\frac{4\gamma }{D}\hfill \end{array}$$

(16)

The friction factor f was computed with the following Equation [37]:

$$f=\left\{\begin{array}{cc}\hfill 64/{\mathrm{Re}}_D\hfill & {\mathrm{Re}}_D\le 3000\hfill \\ \hfill -1.8{log}_{10}\left[{\left({\displaystyle \frac{{\epsilon }_r/D}{3.7}}\right)}^{1.11}+{\displaystyle \frac{6.9}{{\mathrm{Re}}_D}}\right]\hfill & {\mathrm{Re}}_D>3000\hfill \end{array}\right.$$

(17)

where ${\mathrm{Re}}_D$ is the Reynolds number, and ${\epsilon }_r$ is the roughness. The Nusselt number ${\mathrm{Nu}}_D$ was calculated with the following Gnielinski correlation [38].

$${\mathrm{Nu}}_D=\left\{\begin{array}{cc}\hfill 4.36\hfill & {\mathrm{Re}}_D\le 3000\hfill \\ {\displaystyle \frac{(f/8)({\mathrm{Re}}_D-1000)\mathrm{Pr}}{1+12.7\sqrt{f/8}\left({\mathrm{Pr}}^{2/3}-1\right)}}\hfill & {\mathrm{Re}}_D>3000\hfill \end{array}\right.$$

(18)

where $\mathrm{Pr}$ is the Prandtl number. The heat transfer coefficient in Equations (14) and (15) was estimated from the Nusselt number calculated in Equation (18).

2.1.5. Boundary Conditions

The on-gas and on-blast periods times were 80 and 100 min, respectively. In the on-gas period, the fuel gas was supplied with a flow rate of $13.0\times {10}^{-4}$ Nm³/h at 423 K, and its composition is represented in

Table 1. Mass fraction of the fuel gas.

Component	CO2	O2	CO	H2	H2O	N2
Composition (%)	18.1	0.1	24.5	0.3	0.2	56.8

. The air was supplied with a flow rate of $9.2\times {10}^{-4}$ Nm³/h at 383 K. The excess air ratio for the combustion was 1.18. The flue gas outlet was set as a pressure boundary with a gauge pressure of 0 Pa. In the on-blast period, the air was supplied with a flow rate of $36.0\times {10}^{-4}$ Nm³/h at 473 K. The hot blast outlet was set as a pressure boundary with a gauge pressure of 0 Pa. All surface walls were set as nonslip and adiabatic walls. The material properties of the checker bricks that are represented in

Table 2. Material properties of the checker bricks.

Checker Type	Density (kg/m3)	Thermal Conductivity (W/m·K)			Specific Heat Capacity (J/kg·K)
		Temperature (K)			Temperature (K)
		300	1000	2000	300	500	1000	1500	2000
Alumina	2.57	1.70	1.78	1.96	823	994	1171	1249	1284
Fireclay	2.29	1.36	1.48	1.55	799	965	1138	1214	1248

were calculated based on the rule of mixtures and the properties of each component reported by Carnigila and Barna [39].

2.2. Thermal Boundary Conditions of the Structure

During the operation of the hot blast stove, there is heat transfer on the interface between the surface of the structure and the fluids flowing. The types of this heat transfer are divided into convection and radiation. The heat transfer coefficients for both types of heat transfers were calculated from the CFD analysis result. They were provided to the FEA mesh using the user subroutine FILM in Abaqus.

2.2.1. Convective Heat Transfer

The convective heat transfer on the hot blast stove structure was considered with the analytic adiabatic Y-plus method [21]. The method is applicable to CFD calculations based on adiabatic wall boundary conditions. It is possible to analyse the temperature of the structure with conjugate heat transfer; however, it requires enormous computational efforts since numerous cycles have to be calculated to obtain a periodic steady-state solution. In this study, therefore, the temperature of the structure was calculated in an uncoupled scheme to reduce the computational cost with reasonable precision loss [40].

The calculated variables from the CFD were saved into the cell centre, as depicted in Figure 3. The saved variables were the material properties (density, viscosity, and specific heat), the turbulent kinetic energy of the cell centre, the magnitude of the velocity and the fluid temperature. In the standard wall function, the dimensionless cell-wall distance $y^{*}$ and dimensionless velocity $u^{*}$ are calculated as follows, respectively [41].

$$\begin{array}{cc}\hfill y^{*}& =\frac{\rho C_{\mu }^{1/4}{k}_P^{1/2}{y}_P}{\mu }\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill u^{*}& =\left\{\begin{array}{cc}{\displaystyle \frac{1}{\kappa }}ln\left(E{y}^{*}\right)\hfill & y^{*}\ge 11.225\hfill \\ y^{*}\hfill & y^{*}<11.225\hfill \end{array}\right.\hfill \end{array}$$

(20)

In Equations (19) and (20), the subscript P denotes the cell centre, $y_P$ is the distance between the cell centre and the wall, $\kappa$ is the von Kármán constant that is defined as 0.4187, and E is the empirical constant that is defined as 9.793. From the dimensionless velocity $u^{*}$, wall shear stress ${\tau }_w$ and friction velocity $u_{\tau }$ are calculated by the equations as follows:

$$\begin{array}{cc}\hfill {\tau }_w& =\frac{\rho C_{\mu }^{1/4}{k}_P^{1/2}{u}_P}{u^{*}}\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill u_{\tau }& =\sqrt{{\tau }_w/\rho }\hfill \end{array}$$

(22)

where $u_P$ is the magnitude of the velocity of the cell centre, which is perpendicular to the normal direction of the wall. From Kader [42], the temperature profile $T^{+}$ in the fully turbulent boundary layer is expressed as a function of the Prandtl number and dimensionless cell-wall distance of the fully turbulent boundary layer $y^{+}$, which is defined as Equation (23).

$$y^{+}=\frac{\rho u_{\tau }{y}_P}{\mu }$$

(23)

The coefficients $\beta$ and $\mathsf{\Gamma }$ that are required to calculate $T^{+}$ are defined as Equation (24) and Equation (25), respectively [42].

$$\begin{array}{c}\beta ={(3.85{\mathrm{Pr}}^{1/3}-1.3)}^2+2.12ln\mathrm{Pr}\end{array}$$

(24)

$$\begin{array}{c}\mathsf{\Gamma }=\frac{{10}^{-2}{\left(\mathrm{Pr}{y}^{+}\right)}^4}{1+5{\mathrm{Pr}}^3{y}^{+}}\end{array}$$

(25)

$T^{+}$ is calculated as Equation (26) [42].

$$T^{+}\left(y^{+},\mathrm{Pr}\right)=\mathrm{Pr}{y}^{+}+e^{-\mathsf{\Gamma }}+\left[2.12ln\left(y^{+}\right)+\beta \right]e^{-1/\mathsf{\Gamma }}$$

(26)

The convective heat transfer coefficient $h_{\mathrm{conv}}$ is defined as Equation (27), and the convective heat flux ${\dot{q}}_{\mathrm{conv}}^{″}$ is calculated as Equation (28), where $T_{nw}$ is the fluid temperature of the wall adjacent cell, and $T_w$ is the wall temperature [40].

$$\begin{array}{c}{h}_{\mathrm{conv}}=\frac{\rho c_p{u}_{\tau }}{T^{+}}\end{array}$$

(27)

$$\begin{array}{c}{\dot{q}}_{\mathrm{conv}}^{″}=h_{\mathrm{conv}}(T_{nw}-T_w)\end{array}$$

(28)

2.2.2. Radiation Heat Transfer

The heat flux between the surface of the structure and the fluid due to radiation ${\dot{q}}_{\mathrm{rad}}^{″}$ is calculated by Equation (29) as follows:

$${\dot{q}}_{\mathrm{rad}}^{″}=\frac{\sigma \left(\epsilon T_{nw}^4-\alpha T_w^4\right)}{T_{nw}-T_w}$$

(29)

where $\alpha$ is the absorptivity of the grey gas. According to Kirchhoff’s law [43], the emissivity and absorptivity are equal for the opaque surface. Assuming that the surfaces of the refractory linings are opaque, the heat transfer coefficient due to radiation $h_{\mathrm{rad}}$ can be written as Equation (30).

$$h_{\mathrm{rad}}=\sigma \epsilon \frac{T_{nw}^4-T_w^4}{T_{nw}-T_w}$$

(30)

The emissivity of H₂O-CO₂mix gas is calculated using the weighted sum of grey gases model (WSGGM) proposed by Bordbar et al. [44]. The total heat flux ${\dot{q}}^{″}$ on the wall is calculated as Equation (31).

$${\dot{q}}^{″}=(h_{\mathrm{conv}}+h_{\mathrm{rad}})(T_{nw}-T_w)$$

(31)

In the user subroutine FILM, the total heat transfer coefficient $h_{\mathrm{conv}}+h_{\mathrm{rad}}$ and $T_{\mathrm{nw}}$ were provided as the heat transfer coefficient and the sink temperature, respectively.

2.3. Mapping of the Thermal Boundary Conditions

When the nodal positions of meshes used in CFD and FEA are different from each other, the CFD analysis results should be appropriately interpolated and mapped to the FEA mesh. As shown in Figure 4, consider the nodal position s of the mesh to be mapped and the three wall-adjacent cell centres ${\mathbf{x}}_1$, ${\mathbf{x}}_2$, and ${\mathbf{x}}_3$ from the CFD mesh, which are closest from s. The selected cell centres form a plane. Let p is the projection of s onto the plane. All the coordinates of the points are represented in the global coordinate system (${\mathbf{e}}_1$, ${\mathbf{e}}_2$, ${\mathbf{e}}_3$). A new coordinate system (${\mathbf{e}}_1^{\prime }$, ${\mathbf{e}}_2^{\prime }$, ${\mathbf{e}}_3^{\prime }$) can be defined as Equation (32) using ${\mathbf{e}}_1^{\prime }$, which has the same direction as ${\mathbf{x}}_2-{\mathbf{x}}_1$ and the normal of plane ${\mathbf{e}}_3^{\prime }$.

$$\begin{array}{cc}\hfill {\mathbf{e}}_1^{\prime }& =({\mathbf{x}}_2-{\mathbf{x}}_1)/\left|{\mathbf{x}}_2-{\mathbf{x}}_1\right|\hfill \\ \hfill {\mathbf{e}}_3^{\prime }& ={\mathbf{e}}_1^{\prime }\times ({\mathbf{x}}_3-{\mathbf{x}}_1)/\left|{\mathbf{e}}_1^{\prime }\times ({\mathbf{x}}_3-{\mathbf{x}}_1)\right|\hfill \\ \hfill {\mathbf{e}}_2^{\prime }& ={\mathbf{e}}_3^{\prime }\times {\mathbf{e}}_1^{\prime }\hfill \end{array}$$

(32)

In the new coordinate system, $\mathbf{p}-\mathbf{s}$ has the same direction as ${\mathbf{e}}_3^{\prime }$. Translating all the points so that ${\mathbf{x}}_1$ becomes the origin and transforming the coordinate system to the new one, the third-axis coordinates of the cell centres become 0. Setting ${\mathbf{p}}^{\prime }={\left(p_1^{\prime },p_2^{\prime },0\right)}^T={\left(s_1^{\prime },s_2^{\prime },0\right)}^T$, the transformed coordinates of ${\mathbf{x}}_i$ as ${\mathbf{x}}_i^{\prime }={\left(x_i^{\prime },y_i^{\prime },0\right)}^T$, and a scalar to interpolate in i-th cell centre as ${\varphi }_i$, the interpolated scalar value $\varphi$ on p can be calculated using the interpolation function of a linear triangular element in finite element analysis as Equation (33).

$$\begin{array}{c}\hfill \varphi \left(p_1^{\prime },p_2^{\prime }\right)=\frac{1}{J}\left[\left\{x_2^{\prime }{y}_3^{\prime }-x_3^{\prime }{y}_2^{\prime }+(y_2^{\prime }-y_3^{\prime })p_1^{\prime }+(x_3^{\prime }-x_2^{\prime })p_2^{\prime }\right\}{\varphi }_1\right.\\ \hfill +\left\{x_3^{\prime }{y}_1^{\prime }-x_1^{\prime }{y}_3^{\prime }+(y_3^{\prime }-y_1^{\prime })p_1^{\prime }+(x_1^{\prime }-x_3^{\prime })p_2^{\prime }\right\}{\varphi }_2\\ \hfill \left.+\left\{x_1^{\prime }{y}_2^{\prime }-x_2^{\prime }{y}_1^{\prime }+(y_1^{\prime }-y_2^{\prime })p_1^{\prime }+(x_2^{\prime }-x_1^{\prime })p_2^{\prime }\right\}{\varphi }_3\right]\end{array}$$

(33)

In Equation (33), $J=x_2^{\prime }{y}_3^{\prime }-x_3^{\prime }{y}_2^{\prime }+x_3^{\prime }{y}_1^{\prime }-x_1^{\prime }{y}_3^{\prime }+x_1^{\prime }{y}_2^{\prime }-x_2^{\prime }{y}_1^{\prime }$. In this study, the interpolated values at p were used to give the thermal boundary condition on s.

2.4. Material Properties of the Refractory Linings and Shell

The types and applied areas of the refractory linings are described in Figure 5. The apparent density, thermal conductivity, and thermal expansion coefficient of each refractory are given in

Table 3. Material properties of the refractory linings.

Material	Density (g/cm3)	Thermal Conductivity (W/m·K)		Thermal Expansion Coefficient (${10}^{-5}$ K${}^{-1}$)
		at 500 °C	at 1000 °C
Refractory A	2.56	1.77	1.79	64
Refractory B	2.24	1.75	1.79	62
Refractory C	2.11	1.75	1.78	59
Refractory D	0.78	0.24	0.39	55

. The specific heat capacity of the alumina checker brick was used for Refractories A, B and C. For Refractory D, the specific heat capacity of fireclay checker brick was applied. The elastic modulus and Poisson’s ratio of the refractories were set as 4 GPa and 0.15 [45], respectively. The material properties of the refractories were found in the sources [39,45,46] and the catalogues from Chosun Refractories. The shell material was AISI 420 stainless steel with a thickness of 30 mm. Its material properties were obtained from Simufact MatILDa.

2.5. Computational Conditions

In CFD analysis, cyclic analyses were conducted to obtain a periodic steady-state solution for the operation of the hot blast stove. In the cyclic analyses, the state at the end of the on-blast period analysis was used as the initial condition in the on-gas period. The time step was set as 5 s. The velocity and energy were set to be converged if their absolute residuals became less than ${10}^{-3}$ and ${10}^{-5}$, respectively. The periodic steady state was considered to be achieved when the temperature between the cycles became less than 1%.

A grid dependence test was conducted to determine the adequate mesh size for the CFD analysis. Temperature and velocity magnitudes in the middle of the on-gas period were compared with different mesh sizes. Temperature and velocity magnitudes were extracted from the centre of the combustion chamber, as shown in Figure 6a where $x^{*}$ denotes the normalised height. The grid dependence was tested for the mesh sizes of 150 mm, 100 mm, 50 mm and 25 mm. Figure 6b,c show the temperature and velocity profiles of the tested mesh sizes. The temperature and velocity magnitude profiles with the 50 mm mesh size had 0.55% and 0.35% average differences from the 25 mm mesh size, respectively. Considering the accuracy and computational burden, the 25 mm mesh size was selected to generate the mesh for the CFD analysis. The mesh for the CFD analysis was generated with 2,874,424 cells, as shown in Figure 7a.

The finite element analyses were conducted after the CFD analysis achieved its periodic steady-state solution. The thermal boundary conditions were calculated based on the periodic steady-state solution. Structural temperature analyses were conducted to obtain the periodic steady-state solution. The same criterion is used in the structural temperature analysis to decide whether the solution reaches the periodic steady-state solution. The generated FEA mesh has 737,865 elements for the refractory linings and 368,237 elements for the shell, as shown in Figure 7b.

3. Results and Discussion

The velocity profile of the fluid and the temperature profiles of the fluid and checker bricks in the middle of the on-gas period ($t=40$ min) are presented in Figure 8. The maximum temperature reached 1841 K in the on-gas period. During the on-gas period, the velocity profile in the combustion chamber increased to 38.8 m/s at maximum and subsequently decreased to 22.8 m/s at the entrance of the dome due to the cross-sectional area change. A vortex was generated in the centre of the dome while the main flow passed through the outer region of the dome. The vortex caused the main flow to occur through the outer region of the checker chamber. Meanwhile, the inner section of the checker chamber could not absorb the heat from the flue gas due to a lack of flow, even though the vortex temperature was high enough. The external combustion hot blast stove also shows the vortex, but it is located at the top of the checker chamber dome [12]. The dome combustion type has the vortex between the top of the checker bricks and the dome [7,9]. The cause of the vortex generation is the sudden expansion of the flue gas discharged from the combustion chamber [7]. The major difference between the internal combustion type and the other types was the location of the vortex, which disturbed the flow to the inner section of the checker chamber. As a result, the temperature profile of the checker chamber, which is presented in Figure 8c, had an elevated temperature region in the outer region of the checker chamber, which showed a high temperature difference in the horizontal direction. External and dome combustion types do not have a high temperature difference in the horizontal direction since the generated vortex does not interfere with the fluid flow to the checker bricks.

Figure 9 illustrates the velocity and temperature profiles in the middle of the on-blast ($t=130$ min) period. The maximum temperature reached 1804 K, and the average blast temperature reached 1345 K. The vortex was not generated in the on-blast period, and the velocity distribution was considerably even throughout the domain compared to the on-gas period. This flow regime also able to be seen in the dome and external combustion types [9,12]. The fluid velocity started to increase at the entrance of the combustion chamber, and it was hugely increased at the hot blast outlet due to the reduction in the cross-sectional area. The maximum velocity at the hot blast outlet was 117.0 m/s, and the average velocity was 54.2 m/s. The temperature distribution had the elevated temperature region at the outer region of the hot blast stove. The checker bricks in the inner region of the checker chamber received a lower amount of heat during the on-gas period, leading to the lower temperature distribution in the inner region of the hot blast stove. In addition, the dome and external combustion types have comparably slight temperature differences in the horizontal direction, as the checker bricks received a similar amount of heat in the horizontal direction during the on-gas period.

The temperature distributions of the refractory linings in the middle of the on-gas ($t=40$ min) and on-blast ($t=130$ min) periods are described in Figure 10. The inner surface of the refractory linings had similar temperature distribution to the fluid temperature distribution. The outer side of the checker chamber had a higher temperature, whilst the inner region had a lower temperature during the entire operation. The dome region near the outer side of the checker chamber was continuously exposed to a high temperature, as same as the outer side of the checker chamber. From the standpoint of managing the refractory linings, those regions would be a significant area of interest. Most regions of the combustion chamber and the dome had a moderate temperature in the on-gas period, whereas the temperature increased highly in the on-blast period. A noteworthy phenomenon is that the dome region between the combustion and checker chambers showed insignificant temperature change during the operation, which means the region remained in the moderate temperature range during the operation.

Figure 11 depicts the temperature change in the marked region of the combustion chamber wall through the thickness direction during the operation. The region was exposed to the most fluctuating temperature during the operation. The inner surface of the refractory had a temperature change from 1441 K to 1659 K. However, the temperature change reduced below 10 K from 250 mm from the inner surface and vanished at the outer surface. Similar behaviour also showed in the external combustion hot blast stove [12]. The length scale for thermal diffusion $\Delta L_D$ can be estimated as $\Delta L_D=\sqrt{\kappa \Delta t}$ where $\kappa$ is the thermal diffusivity, and $\Delta t$ is the time scale. Taking Refractory A in Table 3 as an example material to compute the length scale for the thermal diffusion, its thermal diffusivity at 1000 °C is 0.875 mm²/s. Its length scale for thermal diffusion is calculated as 97.2 mm with the time scale of $\Delta t=180$ min, which is the operation time for a cycle. The structure temperature analysis results indicated that the temperature difference vanished from 250 mm, which makes the FEA result considered to be reasonable. The shell temperature, therefore, remained steady during the operation thanks to the low thermal diffusivity of the refractories.

The temperature distribution of the shell is depicted in Figure 12. The checker chamber side exhibited the highest temperature distribution due to continuous exposure to the higher temperature of the outer checker chamber. Upon comparing the temperature distributions of the dome and combustion chamber regions, the dome had a lower temperature since the refractory lining was thicker than the combustion chamber. The partition wall side, located between the combustion and checker chambers, was observed to have a significantly lower temperature distribution than the outer regions. It can be understood that the temperature gradient in the circumferential direction comes from the flow characteristics, and the gradient in the height direction comes from the geometric characteristics. Since the conical region exhibited a higher thermal gradient, it was expected to have higher thermal stress.

The maximum and middle principal thermal stress distributions of the shell are illustrated in Figure 13. It is observed that both principal stresses appeared to be higher along the conical region, which was induced by not only the thermal gradient discussed but the geometric characteristics that can raise the stress concentration. However, the external combustion type was reported to have tensile maximum principal stress, and compressive minimum principal stress in the conical region [12], while the internal combustion type had biaxial tensile stress. The conical region in the external combustion hot blast stove has a relatively uniform temperature distribution. Since the thicker part has lower temperature distribution than the thinner part of the conical region, the thinner part tends to expand more than the thicker part. As a result, the conical region of the external type has compressive minimum principal stress, whilst the internal combustion type showed uneven temperature distribution, which generated the tensile minimum principal stress on the conical region of the internal type. An intermediate level of the maximum principal stress was observed on the partition wall side compared to the combustion and checker chamber sides. The stress on the partition wall side was induced due to the temperature gradient in the circumferential direction. The combustion, checker chambers, and partition wall sides had a similar level of middle principal stress since the thermal gradient in the height direction did not differ when compared to the circumferential direction. The maximum principal stress of the dome appeared to be higher in the region near the partition wall side, where a higher thermal gradient was generated.

Figure 14 describes the maximum and middle principal thermal stresses of the conical region along the circumferential direction. The highest maximum principal stress occurred at the 65° point, where the low-temperature area of the conical region was as discussed previously, and its magnitude was 300.6 MPa. The combustion chamber side had higher maximum principal stress distribution than the checker chamber side. A higher middle principal stress appeared in the angle region between 30° and 120°. The highest middle principal stress occurred with 192.0 MPa at the 30° and 105° points. However, the external combustion type has the maximum principal magnitude of 53 MPa [12]. Even though the region which has the highest maximum principal stress is the same in both types, the internal type had about three times higher maximum principal stress. It can be understood that temperature variation is a dominant factor for the generation of thermal stress in hot blast stove applications.

4. Conclusions

In this study, a numerical model was developed to analyse the flow phenomenon, temperature, and thermal stress distribution of a hot blast stove with an internal combustion chamber. The thermal boundary conditions of the hot blast stove structure were calculated using the analytic adiabatic Y-plus method and mapped to the FEA mesh through the developed scheme. The conclusions based on the numerical result can be summarised as follows:

• In the on-gas period, a vortex was generated above the partition wall while the main flow passed through the outer region of the dome. The inner region of the checker chamber could not absorb the heat from the flue gas due to a lack of flow, even though the vortex temperature was high enough;
• The average temperature of the hot blast reached 1345 K in the on-blast period. There was no vortex generation, and the velocity distribution was considerably even compared to the on-gas period. The temperature distribution had the elevated temperature region at the outer region of the hot blast stove;
• The refractory linings in the outer side of the checker chamber and the dome region near the outer side of the checker chamber were continuously exposed to a high temperature. These regions would be a significant area of interest in managing the deterioration of the refractory linings;
• Although the inner surface of the refractory lining temperature changed from 1441 K to 1659 K, the temperature change from 250 mm location into the thickness direction significantly decreased. The shell temperature did not change during the operation thanks to the low thermal diffusivity of the refractories;
• The conical region of the shell was found to have high maximum and middle principal thermal stresses. The maximum and middle principal stresses were calculated as 300.6 MPa and 192.0 MPa, respectively. The thermal stress is induced due to the high-temperature gradient near the conical region, and the geometric characteristic of the conical region made the thermal stress much higher. The thermal stress analysis result indicates that the conical region would be a significant area of interest.

However, there are some limitations in the presented numerical model, which could be suggested in future works. The mechanical stress generated due to operating pressure and gravity was not considered in the model. The CFD model was built with the adiabatic wall boundary conditions to reduce the computational cost. Even if the temperature was calculated with reasonable precision loss, the conjugate heat transfer approach would give a more accurate result. Another limitation is that the degradation of the refractory linings was not considered in the numerical model, where the deterioration worsens the thermal property of the refractories.