Correlations of Heat Transfer and Fluid Flow Data for Lattice Brick Settings in Tunnel Kilns

Abstract

This paper proposes correlation equations of heat transfer and pressure drop for the design and operation of tunnel kilns loaded with lattice brick settings of different geometrical parameters in the form of Nusselt number and friction factor. The developed correlation equations considered parameters that were not investigated in previous studies, such as the relative roughness of the bricks and the stack channels, and they also extended the Reynolds numbers to a practical range that was not covered before in a simple, practical form. The correlation equations are valid for Reynolds number between 125 and 10,200, Prandtl number between 0.68 and 0.73, brick’s relative roughness between 0.23 and 0.93, voidage fraction between 0.48 and 0.653, and the geometrical parameters of the tested lattice brick settings. The achieved correlations of the Nusselt number and the friction factor are well compared with the available correlations in the literature in their valid range of parameters. It is found that Nusselt numbers and the friction factors for low-density are higher than those of high-density settings for all considered parameters except the voidage fraction. The effect of the considered parameters confirms that Nusselt numbers increase and the friction factors decrease substantially with the Reynolds number and slightly with the Prandtl number. At a constant Reynolds number, both the Nusselt number and the friction factor increase as the brick’s relative roughness is increased. Moreover, as the stack channel spacing is increased, the Nusselt number decreases, and the friction factor increases. The voidage fraction of the setting has a monotonic effect on both Nusselt numbers and friction factors. Nusselt numbers for high-density are higher than those for low-density settings as the voidage fraction varies.

1. Introduction

Ceramic materials manufacturing is a huge industry worldwide that consumes a great amount of energy. Bricks, wall and floor tiles, roof tiles, sanitary ware, vitrified clay products, and table-and-ornamental wares are examples of ceramic products that are made of clay minerals. Ceramic products can be dried and fired using many types of furnaces, including tunnel kilns which can be divided into three main temperature zones; preheating, firing, and cooling zones. Tunnel kilns can be considered counter-flow heat exchangers because the raw ceramic materials enter the tunnel in the preheating zone where the combustion gases exit. A great amount of fuel is consumed in the firing zone to achieve high temperatures, which is required for brick sintering [1,2].

Depending on the type of manufactured bricks, the firing temperature could vary from 800 to 1800 °C. For example, common bricks require a firing temperature in the range of 900–1050 °C, while basic bricks require a firing temperature in the range of 1400–1800 °C [3]. Such high temperatures make the brick manufacturing industry an energy-intensive industry. With the specific energy consumption (SEC) being defined as the amount of energy used to produce 1 kg of fired bricks, the actual brick SEC ranges from 2.31 to 3.51 MJ/kg, depending on the fuel type and kiln efficiency [4]. Other researchers reported different rates of SEC, between 2.04 and 3.51 MJ/kg [5] and 3.47 MJ/kg [6]. Besides, SEC could be higher for ineffective kilns. As the actual specific energy consumption in tunnel kilns is more than double the theoretical one, many studies attempt to understand and analyze its thermal performance.

Many studies have been performed in the literature to investigate energy consumption and thermal management in tunnel kilns. Riedel [7] found that convection is the key factor in analyzing heat transfer and getting optimum results for tunnel kilns. The first recognized convection heat transfer data in the literature was developed by Abbakamow [8], which was criticized by Dugwell and Oakley [9] as it predicts half the Nu predicted using their correlation.

Then, Dugwell and Oakley [9] examined the heat transfer in tunnel kilns experimentally and proposed their correlation. They simulated airflow through two different ware settings that include three and five blades and estimated the convective heat transfer coefficients for the chrome magnesite block. However, the brick columns in their study were treated as one solid column without stack channels. Based on their approach, the airflow in the column channels was not investigated, and this could be considered a drawback in their correlation. In addition, Dugwell and Oakley [9] developed different correlations for the different faces of the blade (top, front, back, and side faces). Thus, it could be assumed as a local correlation for column faces, not the brick settings.

Abou-Ziyan [10] conducted an experimental investigation to study the thermal performance of six different brick setting arrangements. The setting pattern was found to be a significant factor that affects the pressure drop and convective heat transfer coefficient greatly. Carvalho and Nogueira [11] utilized a modeling tool to optimize several thermal equipments such as cement tunnel kilns, ceramic tunnel kilns, baking ovens, and glass-melting furnaces for different heat transfer processes. Based on their study, a substantial reduction in energy consumption and pollution has been accomplished. This showed that numerical tools could be considered effective tools in optimization procedures for kiln design and operation.

The lattice brick setting was found to affect the heat transfer and brick production in tunnel kilns [10,12,13,14,15,16,17]. The various tested settings in those studies have different spacing between the bricks in a single column, which is named stack channels in the present study. Thus, stack channels play a key role in the production of tunnel kilns and the quality of the bricks. Furthermore, a 3D-CFD model was used to investigate the airflow uniformity within the brick setting and the convective heat transfer coefficients for a specific lattice setting for a different column, wall, ceiling, and extension channel dimensions [18,19]. Optimum sizes of column, wall, ceiling, and extension channels were suggested among different tested dimensions that produce the best flow uniformity and heat transfer between bricks and flowing air. Dmytrochenkova and Tadlya [20] carried out 2D numerical simulations to investigate the airflow distributions in the brick-setting channels placed in tunnel kilns. They found that the reverse airflow in the cooling zone can be minimized, and the heat loss can be reduced by reducing the arch height. However, the pressure drop was found to be increased. Almesri et al. [21] and Alrahmani et al. [22] performed 3D numerical studies to investigate the effect of surface roughness on fluid flow and heat transfer in tunnel kilns. They found that the surface roughness plays a crucial role in enhancing the heat transfer in tunnel kilns and hence the quality of brick production. Alrahmani et al. [23] studied the effect of brick setting densities on the airflow pattern and heat transfer numerically by comparing two different types of lattice brick settings of different densities. They found that the setting of fewer bricks (low density) enhances the movement of airflow within the brick setting, especially in the stack channels when compared with the high-density setting and hence the heating of the setting bricks was augmented. Furthermore, the pressure drop and pumping power were reported to be less in the case of low-density settings, which consume less energy.

Many studies have been conducted to propose empirical correlation equations of Nusselt number and friction factor, which are beneficial in designing tunnel kilns and provide the designers with accurate predictions of heat transfer and fluid flow characteristics. For example, Abou-Ziyan [10] obtained correlation equations to predict the pressure drop and heat transfer for refractory materials in the form of friction factor and Nusselt number as functions of Reynolds number, voidage fraction (ε) (the ratio of the free volume of the wind tunnel to the effective tunnel volume), column channel spacing to brick thickness (CC/a) and column channel spacing to brick length (CC/c). The correlations are as follows:

$$f={\left(\frac{CC}{a}\right)}^{-5.18}{\left(\frac{CC}{c}\right)}^{5.83}{\epsilon }^{1.68}\left\{\frac{7010}{R{e}^{0.2}}-300{\left(\frac{CC}{a}.\frac{CC}{c}\right)}^{0.53}\right\}$$

(1)

$$Nu=0.088 R{e}^{0.625}{\left(\frac{CC}{c}\right)}^{-2.30}{\left(\epsilon \frac{CC}{a}\right)}^{1.837}$$

(2)

The correlation was validated for Re range from 6000 to 33,000, voidage fraction ranges from 0.43 to 0.57, (CC/c) from 0.44 to 1.17, and (CC/a) from 1.62 to 4.25.

Vogt and Beckmann [24] obtained a correlation equation to calculate the coefficient of friction in tunnel kilns. The correlation helped in understanding the amount of pressure loss in extension channels of brick settings placed in a tunnel. The correlation is dependent on the ratio of channel length to the equivalent diameter, the hydraulic diameter of the rectangular tunnel (G₁), the ratio of the extension length to the equivalent diameter (G₂), and the ratio of channel area to the extension area (G₃).

$$f=-0.4+0.48·\left(1-e^{-4G_2}\right)+\frac{300}{G_2^{5.0}·\left(e^{\frac{20}{G_2}-33·G_2^{0.62}·e^{-G_2^{1.35}}·G_3^{0.457-5.9·G_2^4·e^{-3·G_2}}}-1\right)} +\left(\frac{64}{Re}\right)·\left(0.1+\left(1-e^{-\left(0.028+\left(1-e^{-10000·G_3^{9.0}}\right)·2.4\right)·G_1^4}\right)·4.4\right)$$

(3)

The equation is valid over the following ranges: 30 < Re < 19,000, 0 < G₁ < 12, 0.08 < G₂ < 4.9, and 0.1 < G₃ < 0.67.

Vogt and Beckmann [25] continued their study on understanding the pressure loss, flow distribution, and convective heat transfer in brick settings in tunnel kilns. They proposed generalized correlation equations for the Nusselt number by utilizing the measurements of Dugwell and Oakley [9]. They improved the correlation equation of Abbakumow [8], which was found to be inadequate as it does not include the effect of the channel flow in the extension properly. Vogt and Beckmann [25] utilized in their correlation equations the three non-dimensional similarity indices (G₁, G₂, and G₃), which were introduced in their previous studies [24]. They found that G₁ and G₂ equivalently affect the convective heat transfer in the channel. While G₁ revealed a more significant effect on the convective heat transfer in the extension compared to the effect of G₂ and G₃. For the channels, they used Abbakumow’s correlation equation for the Nusselt number [8]:

$$Nu=\mathsf{\Phi }·0.018·R{e}^{0.8}.$$

(4)

Furthermore, from the measurements of Dugwell and Oakley [9], they proposed the parameter (Φ) as follows:

$$\mathsf{\Phi }=1.64+0.71·\left(e^{-0.035G_1^{1.8}}-1\right)+0.17\left(1-e^{-1.1G_1}\right)+1.05\left(1-e^{-1.5\times {10}^7{G}_3^{18.0}}\right)+\left(2.36+2.08\left(e^{-0.46G_1}\right)-1\right)\left(1-e^{-\left(0.23+0.5\left(1-e^{-0.0007G_1^{5.0}}\right)\right)G_2^2}\right)-\left(-0.31-0.225\left(1-e^{-0.52G_1}\right)\right)+\left(1.03-0.94\left(1-e^{-0.52G_1}\right)\left(2\left(1-e^{-0.35·G_2^{1.5}}\right)\right)\left(1-e^{\left(-\left(10-4\left(1-e^{-0.9G_2}\right)\right)-1.795e^{-\frac{{\left(G_2-0.85\right)}^{2.0}}{0.08}}\right)\times {10}^{-5}Re}\right)\right)$$

(5)

The equation extends over the following fields: 6100 < Re < 50,700, 0.1 < G₁ < 5.2, 0.36 < G₂ < ∞, and 0.3 < G₃ < 0.65.

Refaey et al. [16] enhanced the rate of heat transfer, in the cooling zone, by adding guide vanes in the tunnel kiln for ten different brick settings. In addition, a Nusselt number empirical correlation in terms of Reynolds number (Re), brick setting dimensionless groups (CC/c), (ε CC/a), and the guide vane’s angle of attack (θ) was obtained. Equation (6) is valid for Re between 11,867 and 25,821.

$$Nu=1.617 R{e}^{0.42}{\left(\frac{CC}{c}\right)}^{-0.147}{\left(\mathsf{\epsilon }\frac{CC}{a}\right)}^{0.224}{\left(\frac{\theta }{180}\right)}^{-0.39}$$

(6)

where CC represents the spacing between columns, c is the brick length, a is the brick thickness, and ε is the voidage fraction. The correlation satisfied the collected experimental data within ±15% maximum deviation.

Refaey et al. [17] suggested another approach to enhance heat transfer in tunnel kilns by using two different guide vane types, namely, U-shape guide vanes and side wall guide vanes. They obtained experimental correlations for vane types for the average Nusselt number in terms of Reynolds number (Re), brick setting dimensionless group (CC/c), ($\mathsf{\epsilon }$ CC/a), and the guide vane’s angle of attack ($\theta$). Their correlations are valid for Re between 13,609 and 27,634.

$$Nu=1.141 R{e}^{0.427}{\left(\frac{CC}{c}\right)}^{-0.313}{\left(\mathsf{\epsilon }\frac{CC}{a}\right)}^{0.255}{\left(\frac{\theta }{180}\right)}^{-0.639}(\mathrm{for}\mathrm{U}-\mathrm{shape}\mathrm{guide}\mathrm{vanes})$$

(7)

$$Nu=2.362 R{e}^{0.358}{\left(\frac{CC}{c}\right)}^{-0.259}{\left(\mathsf{\epsilon }\frac{CC}{a}\right)}^{0.213}{\left(\frac{\theta }{180}\right)}^{-0.346}\left(\mathrm{for}\mathrm{side}\mathrm{wall}\mathrm{guide}\mathrm{vanes}\right)$$

(8)

The above literature survey revealed that although the recent empirical correlations of Nusselt number and friction factor are correlated with parameters such as Reynolds number, voidage fraction, column channels, and extension channels. However, those correlations did not consider some important parameters, such as the effect of brick surface roughness and the stack channels. Moreover, the correlations are applicable only for high Reynolds numbers. Therefore, the present study was performed using CFD to develop correlation equations for the Nusselt number and friction factor with parameters such as Reynolds number, Prandtl number, voidage fraction, relative roughness of bricks, the ratio of column channel spacing to brick length $\left(\frac{CC}{c}\right)$, the ratio of extension channel spacing to brick height $\left(\frac{EC}{b}\right)$ and the ratio of stack channel spacing to brick thickness $\left(\frac{SC}{a}\right)$ for two different lattice setting density, which is classified by their voidage fraction, namely, low density-setting (LDS) and high-density setting (HDS). The LDS is composed of 512 bricks and has a voidage fraction of 0.653, while the HDS consists of 768 bricks with a voidage fraction of 0.48. The proposed correlation equations in the present study are also applicable to low Reynolds numbers. The developed correlations provide reliable design and operation data for tunnel kiln developers and practitioners.

2. Numerical Simulation

This section presents the numerical modeling for the simulation of airflow and heat transfer in a tunnel kiln loaded with two different lattice brick settings: LDS and HDS, in order to develop correlation equations for both Nusselt number (Nu) and the friction factor (f) with Reynolds number (Re), Prandtl number (Pr), relative roughness of bricks (RR), voidage fraction ($\epsilon$), $\left(\frac{CC}{c}\right)$, $\left(\frac{EC}{b}\right)$ and $\left(\frac{SC}{a}\right)$ where $CC$ is the column channel spacing, EC is the extension channel spacing, SC is the stack channel spacing, a is the brick thickness, b is the brick height, and c is the brick length (the geometrical parameters a, b and c are shown in Figure 1, while the column channels, extension channels, and stack channels are illustrated in Figure 2). The commercial Computational Fluid Dynamics (CFD) code by ANSYS Fluent software version 18.2 is used to construct the mathematical model in which the finite volume method is implemented to solve the governing equations (i.e., the mass, momentum, and energy conservation) numerically. The computational domain was discretized and solved using an unstructured grid with refined cells near solid surfaces, i.e., the bricks and the tunnel walls. The second-order upwind scheme is used for the pressure, momentum, energy, turbulent kinetic energy, and specific dissipation rate. The pressure and velocity fields are combined using the “coupled” algorithm. The prescribed convergence criterion for the equations of continuity and momentum is set to be ${10}^{-7}$, while for the energy equation is ${10}^{-8}$.

The governing equations, which are solved in the computational cells, are the conservation of mass, conservation of momentum, and conservation of energy. These governing equations are represented by Equations (9)–(11), respectively.

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}\left(\rho u_i\right)=0$$

(9)

where $\rho$ is the density of air, $t$ is the time, $u_i$ is the velocity vector components (u, v, and w), and $x_i$ is the Cartesian coordinate axis (x, y, and z).

$$\frac{\partial }{\partial t}\left(\rho u_i\right)+\frac{\partial }{\partial x_j}\left(\rho u_i{u}_j\right)=-\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)\right]+\rho g_i$$

(10)

where p represents the pressure, $\mu$ is the dynamic viscosity, and $\rho g_i$ is the body force in the x, y, and z directions.

$$\frac{\partial }{\partial t}\left(\rho T\right)+\frac{\partial }{\partial x_i}\left(\rho u_{i}T\right)=\frac{\partial }{\partial x_i}\left[\frac{k}{c_p}\frac{\partial T}{\partial x_i}\right]+S_T$$

(11)

where T is the temperature, k is the thermal conductivity, $c_p$ is the specific heat, and $S_T$ is the source term. It has to be noted that all properties used in the CFD calculations are a function of temperature. Moreover, the effect of radiation was not considered in the present model.

Time-averaging the above-mentioned governing equations results in steady Reynolds Averaged Navier-Stockes (RANS) equations. The Reynolds stresses, which are derived from the averaging operation, are solved using a turbulence model. Section 2.2 illustrates the ability of the selected turbulence model to predict the airflow pattern across such complex settings as well as the surface temperature of the bricks accurately.

2.1. Computational Domain and Boundary Conditions

In this study, the wind tunnel experiment on a lattice brick set which was performed by Abou-Ziyan [10], is simulated using CFD. Two different lattice brick settings are used in the current study, namely the low and high-density settings. The lattice setting is made of eight rows, and each row consists of four columns. The low-density setting (LDS) and the high-density setting (HDS) are composed of 512 and 768 bricks, respectively. The bricks are made of fire clay, which has a density of 2645 $\mathrm{kg}/{\mathrm{m}}^3$, specific heat of 960 $\mathrm{J}/\mathrm{kg}·\mathrm{K}$ and thermal conductivity of 1.31 $\mathrm{W}/\mathrm{m}·\mathrm{K}$. The dimensions of each brick are 58 × 28 × 16 mm. Some bricks are arranged longitudinally, and others transversely to the direction of the flow. The wind tunnel used in the CFD simulations, as shown in Figure 1, is 3.75 m long, 0.35 m wide, and 0.236 m high. To eliminate the entrance effect, the brick setting is located 1.75 m from the entrance. The arrangement of the brick setting in the kiln tunnel creates five different types of flow channels, namely three column channels (spacing between brick columns), two wall channels (spacing between the wall of the tunnel and the brick column), one ceiling channel (spacing between the ceiling of the tunnel and top layer of the bricks), seven extension channels (spacing between brick rows) and many stack channels (channels between bricks within the columns in the setting), as shown in Figure 2. The dimensions of the ceiling and extension channels for both brick settings are 12 and 21 mm, respectively, while the dimensions of the wall, column, and stack channels are 15.2, 29.2 and 26 mm, respectively, for LDS and 14, 30 and 5 mm, respectively for HDS. These dimensions of different flow channels of the brick setting were found to produce uniform airflow throughout the main tunnel kiln channels [18,19]. Both LDS and HDS are shown in Figure 3.

In order to get accurate CFD results, precise use of the boundary conditions is required. To develop the correlation equations for both the Nusselt number (Nu) and the friction factor (f), thirty simulation runs have been performed for each brick setting (i.e., LDS and HDS) in which the values that are used for the inlet air velocity are ranging from 1 m/s to 9 m/s. Although the used inlet air velocity is the same for both LDS and HDS settings, the Reynolds number is different for both settings due to their different hydraulic diameters and voidage fractions. The values used for the inlet air temperature range from 300 to 2200 K, and the values used for the brick surface roughness range from 1 to 4 mm (relative roughness of bricks (RR): 0.231–0.924). These values of the inlet air velocity, inlet air temperature, and brick surface roughness have been selected using the design of the experiment method. The test conditions of the sixty CFD simulation runs, which are used to develop the correlation equations for the Nusselt number (Nu) and the friction factor (f), are shown in

Table 1. The test conditions used in the CFD simulations for both the LDS and HDS settings.

Test Condition	Inlet Velocity (m/s)	Brick Surface Roughness (mm)	Inlet Temperature (K)
1	1.000	3.85	1535.0
2	1.400	1.45	838.3
3	1.667	1.65	2200.0
4	1.933	2.75	1028.3
5	2.200	2.25	1155.0
6	2.467	2.05	1661.7
7	2.733	3.15	300.0
8	3.000	2.85	1408.3
9	3.267	3.05	1471.7
10	3.533	1.55	1598.3
11	3.800	1.95	521.7
12	4.067	2.65	965.0
13	4.333	1.75	1725.0
14	4.600	3.65	1091.7
15	4.867	1.0	2105.0
16	5.133	3.75	2041.7
17	5.400	1.35	395.0
18	5.667	4.0	1978.3
19	5.933	3.45	1218.3
20	6.200	2.15	711.7
21	6.467	2.45	585.0
22	6.733	3.35	458.3
23	7.000	1.85	775.0
24	7.267	2.35	901.7
25	7.533	1.15	1851.7
26	7.800	2.55	648.3
27	8.067	3.55	1788.3
28	8.333	3.25	1281.7
29	8.600	1.25	1915.0
30	9.000	2.95	1345.0

.

To calculate the convective heat transfer coefficients (CHTCs) in the current study, all longitudinal bricks that are located at all layers of the first and second columns and all transverse bricks that are located at all layers of the third and fourth columns are changed to copper. Those copper bricks are called “instrumented bricks”. The heat rate assigned to the longitudinal and transverse instrumented bricks is 14 W. The values of the density, specific heat and thermal conductivity for the copper bricks which are used in the CFD simulations are 8978 $\mathrm{kg}/{\mathrm{m}}^3$, 381 $\mathrm{J}/\mathrm{kg}·\mathrm{K}$ and 387.6 $\mathrm{W}/\mathrm{m}·\mathrm{K}$, respectively. The tunnel walls are assumed to be adiabatic surfaces. No-slip conditions at the boundaries of the wind tunnel and the bricks surfaces are applied. At the outlet of the wind tunnel, a pressure of zero Pascal is applied.

2.2. Verification and Validation of the CFD Model

To achieve reliable CFD predictions, it is essential to employ a numerical grid with adequate refinement, ensuring that the obtained results are not affected by the grid utilized. Therefore, the wind tunnel meshed using four different grid resolutions with the aim of evaluating grid independence and identifying the most suitable grid that yields a converged solution, accurate results, and a manageable computational timeframe. The pressure drop across the brick setting (from x = 1.7 m to x = 3.0 m), as well as the surface temperature of the longitudinal and transverse bricks, were used to evaluate the predicted values using various grids, as detailed in

Table 2. Grid sensitivity analysis using pressure drop ($\Delta \mathit{P}$) and surface temperature (${\mathit{T}}_{\mathit{s}}$) of the longitudinal and transverse instrumented copper bricks for both the LDS and HDS.

	Computed Values				% Difference
Parameter	Grid 1 (0.72 MC)	Grid 2 (1.44 MC)	Grid 3 (4.78 MC)	Grid 4 (6.64 MC)	Grid 1	Grid 2	Grid 3	Grid 4
$\Delta P$ (pa)	228.97 (330.3)	217.07 (318.5)	213.8 (317.9)	213.58 (313.4)	7.21 (5.39)	1.63 (1.63)	0.10 (1.44)	0 0
Ts (K), long brick	326.04 (322.8)	325.89 (322.5)	325.1 (322.4)	325.06 (322.1)	0.30 (0.22)	0.26 (0.12)	0.012 (0.09)	0 0
Ts (K), trans brick	327.02 (324.5)	326.68 (324.3)	326.5 (323.9)	326.17 (323.8)	0.26 (0.22)	0.16 (0.15)	0.10 (0.03)	0 0

Note: The values of the parameters between brackets are for the high-density setting (HDS).

. Grid sizes ranged from 0.72 Mega cells to 6.64 Mega cells. In general, the disparity between each grid’s predicted value and Grid 4′s value diminished as grid resolution increased. Among the numerical grids examined in this research, Grids 3 and 4 revealed almost grid-independent results compared to the other grids. Since the difference between the predicted values of pressure drop and surface temperatures of the bricks that are produced by Grids 3 and 4 are very small, and Grid 4 requires more computational time, Grid 3 was chosen for this study.

To validate the results of the CFD model, Grid 3 was incorporated with three different turbulence models, namely Standard k-ε, RNG k-ε, and k-ω. The calculated pressure drop and the surface temperature of the longitudinal and transverse instrumented bricks were compared with the experimental data obtained by the Abou-Ziyan study [10]. It is worth mentioning that an inlet air velocity of 6.9 m/s, inlet air temperature of 303 K, and heat generations of 346,165 and 344,231 W/${\mathrm{m}}^3$ for longitudinal and transverse instrumented bricks, respectively, were used in the experimental work by Abou-Ziyan [10] and hence they were used as boundary conditions in the CFD simulations. Figure 4 shows the locations of the longitudinal and transverse copper bricks in the brick setting, which are used for the CFD model validation. The brick settings used in the CFD validation process consisted of only six rows to coincide with the experimental study carried out by Abou-Ziyan [10], and then the brick setting was extended to eight rows to investigate the effect of row numbers on the fluid flow development and the establishment of the heat transfer.

Grid 3 was chosen, as mentioned above, to compare the performance of three different turbulence models (standard k-ε, RNG k-ε, and k-ω) in predicting the characteristics of the airflow across the brick setting and the heat transfer between the flowing air and the bricks’ surfaces for both settings the LDS and HDS.

As shown in

Table 3. The experimental and computed pressure drop ($\Delta \mathit{P}$) and surface temperature (${\mathit{T}}_{\mathit{s}}$) of the longitudinal and transverse instrumented copper bricks for the three tested turbulence models for both the LDS and HDS using Grid 3.

Parameter	Experimental Values	Computed Values			% Difference
		Standard $\mathbf{k}-\mathit{\epsilon }$	$\mathbf{RNG}\mathbf{k}-\mathit{\epsilon }$	$\mathbf{k}-\mathit{\omega }$	$\mathbf{Standard}\mathbf{k}-\mathit{\epsilon }$	RNG $\mathbf{k}-\mathit{\epsilon }$	$\mathbf{k}-\mathit{\omega }$
$\Delta P$ (pa)	210.7 (312)	172.65 (261.6)	140.51 (240.8)	213.80 (317.9)	−18.06 (−16.2)	−33.31 (−22.8)	1.47 (1.89)
Ts (K), long brick	323 (319)	325.10 (322.7)	329.80 (324.7)	325.10 (322.4)	0.65 (1.16)	2.11 (1.79)	0.65 (1.066)
Ts (K), trans brick	325 (322)	324.90 (324.2)	330.71 (326.7)	326.50 (323.9)	−0.031 (0.68)	1.76 (1.46)	0.46 (0.59)

Note: The values of the parameters between brackets are for the high-density setting (HDS).

, the k-$\omega$ turbulence model produced better predictions of pressure drop when compared with the other turbulence models, and the difference between the computed and experimental values of $\Delta P$ of LDS and HDS is only 1.47 and 1.89%, respectively. For the brick surface temperatures, both the Standard k-$\epsilon$ and k-$\omega$ models showed similar performance in predicting the surface temperature of the longitudinal brick for the LDS setting (the percentage difference between computed and measured values is 0.65%), while the Standard k-$\epsilon$ model predicted a slightly better value of the surface temperature of the transverse brick for LDS when compared with k-$\omega$ model (the percentage differences between computed and measured values using Standard k-$\epsilon$ and k-$\omega$ models are 0.031 and 0.46%, respectively). For HDS, k-$\omega$ was found to be the best model for predicting the surface temperature of the bricks (the percentage differences between computed and measured values of the surface temperature of longitudinal and transverse bricks are 1.066 and 0.59%, respectively. The RNG k-$\epsilon$ model showed the worst results among the three tested models. Therefore, the k-$\omega$ model was used to predict the airflow pattern across the brick setting and the bricks’ surface temperatures under various conditions.

3. Data Analysis and Reduction

The convection heat transfer represents over 80% of heat transfer during brick firing in tunnel kilns [25]. The convective heat transfer coefficient between the whole brick row and the flowing air is calculated using Equation (12):

$$h=\frac{\dot{Q}}{A_b\left(T_b-T_a\right)}$$

(12)

where $\dot{Q}$ is the total heat dissipated by the instrumented bricks in the row, $A_b$ is the surface area of all bricks participating in the heat dissipation, $T_b$ is the surface temperature of all instrumented bricks in a certain row, and $T_a$ is the arithmetic average of the air temperature calculated for planes located in the middle of the extension channels before and after each considered row.

The Nusselt number for each brick row is calculated using Equation (13):

$$Nu=\frac{h{D}_h}{k}$$

(13)

where h is the convective heat transfer coefficient of the brick row (Equation (12)), k is the thermal conductivity of the flowing air, which represents the arithmetic average of thermal conductivity that is calculated for planes located in the middle of the extension channels before and after each considered row and $D_h$ is the hydraulic diameter which is calculated using Equation (14):

$$D_h=\frac{4V_f}{A_w}$$

(14)

where $V_f$ is the free volume of the loaded wind tunnel, which is calculated by subtracting the volume of the total bricks ($V_b$) from the effective loaded tunnel volume ($V_t$) (the volume of the wind tunnel that surrounds the brick setting only and without the extension channels),

$$V_f=V_t-V_b$$

(15)

and $A_w$ is the wetted area which is calculated using Equation (16):

$$A_w=A_{b,w}+A_{t,w}$$

(16)

where $A_{b,w}$ is the wetted area of the bricks and $A_{t,w}$ is the wetted area of the loaded tunnel.

The friction factor for each brick row is calculated using Equation (17):

$$f=\frac{\Delta P{D}_h}{\frac{1}{2}\rho L{U}^2}$$

(17)

where $\Delta P$ is the pressure drop across each brick row, $D_h$ is the hydraulic diameter (Equation (14)), $\rho$ is the air density which represents the arithmetic average of air density that is calculated for planes located in the middle of the extension channels before and after each considered row, L is the width of the brick row, and U is the interstitial velocity (the velocity of air within the brick setting which is calculated by dividing the inlet air velocity (V) by the voidage fraction ($\epsilon$) (Equation (18)).

$$\epsilon =\frac{V_f}{V_t}$$

(18)

The Reynolds number for each brick row is calculated using Equation (19):

$$Re=\frac{U{D}_h}{\nu }$$

(19)

where $\nu$ is the kinematic viscosity which represents the arithmetic average of kinematic viscosity that is calculated for planes located in the middle of the extension channels before and after each considered row.

The Prandtl number for each brick row is calculated using Equation (20):

$$Pr=\frac{\mu C_p}{k}$$

(20)

where $C_p$ is the specific heat of air. It is worth mentioning that the values of $\mu$, $C_p$ and k used in Equation (20) represent the arithmetic average that is calculated for planes located in the middle of the extension channels before and after each considered row.

The relative roughness of the bricks RR is calculated using Equation (21):

$$RR=\frac{e_b}{L_c}$$

(21)

where e_b is the roughness of the bricks and $L_c$ is the characteristic length which is calculated by dividing the volume of a single brick (V_br) by its surface area (A):

$$L_c=\frac{V_{br}}{A}$$

(22)

4. Results and Discussion

This section presents the developed Nusselt number and friction factor correlation equations for LDS and HDS. Moreover, general correlation equations for both parameters are developed for both settings, which can serve for the studied parameters of Reynolds number, Prandtl number, relative roughness, and the geometrical parameters of both settings. Such correlations facilitate better design and operation data for tunnel kilns that consider the brick’s relative roughness, which was not considered in the earlier literature correlations of tunnel kilns. Moreover, the parameters of the correlations are discussed to explore their effects on heat transfer and pressure drop in tunnel kilns.

4.1. Correlation of Heat Transfer and Pressure Drop Data

Figure 5 illustrates the variations of Nusselt numbers and friction factors with the row numbers of the bricks. Both Nu and f decrease from row number one to five, where it is stabilized to an approximately constant value. Therefore, the correlation of the Nusselt number and the friction factor is conducted on the average values of rows five to seven. The first few brick rows are influenced by the large turbulence due to the entrance effect. Moreover, row eight is excluded as it sometimes suffers from the exit effect.

The heat transfer and pressure drop data sets for LDS and HDS are simulated for different conditions of interstitial Reynolds number, Re (125–10,200), Prandtl number, Pr (0.68–0.73), and brick relative-roughness, RR (0.23–0.92). Sixty runs are simulated, and the data reduction technique is applied to obtain Re, Nusselt number (Nu), and friction factor (f), as discussed in Section 2. The Nu and f parameters for LDS and HDS are correlated as functions of Re, Pr, and RR using the non-linear least square regression method. In addition, the data for both LDS and HDS are combined and correlated with a single equation for either Nu or f. These general correlations provide reliable data for the efficient design and operation of tunnel kilns to avoid kiln oversize and excessive energy consumption.

4.1.1. Correlations for LDS

Nusselt numbers and friction factors are correlated for LDS, as given in Equations (23) and (24), with R² of 0.9942 and 0.6587, respectively. These equations are valid for wide ranges of Re (185–10,200), Pr (0.68–0.73), and RR (0.23–0.92). Since the variations in Pr are limited, the values of Nu are divided by Pr^0.333, and those of friction factors are divided by Pr^−0.666 before conducting the correlation process. Figure 6a presents the prediction against the simulated Nu, and Figure 6b shows the prediction against simulated f. Most Nu predictions lie within ±5% (Figure 6a), whereas f predictions are within ±10% (Figure 6b). The more accurate predictions of Nu than f are supported by higher R² for Nu than f.

$$Nu=0.039 R{e}^{0.858} P{r}^{0.333} R{R}^{0.26},{\mathrm{R}}^2=0.9942$$

(23)

$$f=0.5672 R{e}^{-0.0545} P{r}^{-0.666} R{R}^{0.36}, R^2=0.6587$$

(24)

In addition to the above-stated correlations, Nu is correlated with the friction factor and Re with a high coefficient of determination (R² = 0.9954), as given in Equation (25). Therefore, Equation (25) provides an accurate evaluation of the friction factors than Equation (24). Figure 6c shows the relation between Nu and $R{e}^{0.831} f^{0.6869}$, where the line has a slope of 0.301. All the data points are distributed around the line representing the relation between Nu and f.

$$Nu=0.07231 R{e}^{0.831} f^{0.689},R^2=0.9954$$

(25)

4.1.2. Correlations for HDS

Similar to the correlations developed for the LDS, Nu and f data sets for the HDS are correlated as given in Equations (26) and (27), with coefficients of determination R² of 0.9953 and 0.8257, respectively. Equations (26) and (27) are valid for Re (125–7600), Pr (0.68–0.73), and RR (0.23–0.92). The friction factors for the HDS are correlated better than those of the LDS. The accurate predictions are indicated by higher R² for the case of HDS (0.8257) than the LDS (0.6587). Figure 6d,e show the distribution of the predicted Nu and f, respectively, against the simulated values. The values of Nu distributed in the envelope of ±5% and those of f are within ±10%, indicating accurate predictions of Nu than f.

$$Nu=0.01 R{e}^{1.0} P{r}^{0.333} R{R}^{0.2526},R^2=0.9953$$

(26)

$$f=0.2616 R{e}^{-0.0673} P{r}^{-0.666} R{R}^{0.363}, R^2=0.8257$$

(27)

Again, Equation (28) presents the correlation between Nu and f and Re for the HDS, with a higher coefficient of determination than Equation (27). Moreover, Figure 6e shows the correlation between Nu and $R{e}^{0.9775} f^{0.8145}$, where the line has a slope of 0.222. Clearly, all the data points lie on the line representing the stated correlation (Figure 6e). More accurate friction factors can be obtained using Equation (28) than Equation (27).

$$Nu=0.04086 R{e}^{0.9775} f^{0.8145},{\mathrm{R}}^2=0.9989$$

(28)

4.1.3. General Correlations for Both Settings

The simulated Nusselt numbers and friction factors for both LDS and HDS are correlated with a single equation for both settings. In order to combine both settings, four geometrical dimensionless groups are introduced to characterize the setting arrangement. Those groups are the kiln voidage fraction (ε), column channel spacing (CC) to brick length (c), extension channel spacing (EC) to brick height (b), and stack channel spacing (SC) to brick thickness (a).

The obtained general correlation equations are valid for Re between 125 and 10,200, Pr from 0.68 to 0.73, relative brick roughness between 0.23 and 0.92, kiln voidage fraction between 0.48 and 0.65, column channel group of about 0.51, extension channel group of 0.75, and stack channel group between 0.31 and 1.31. Nusselt number is well correlated with R² higher than 0.99 as given in Equation (29). Moreover, the friction factor is well correlated with R² of 0.9441, as in Equation (30). Figure 7a,b show the predicted and simulated Nu and f, respectively. Most of the Nu predictions are within ±5%, as shown in Figure 7a, whereas Figure 7b shows that the friction factor predictions are scattered within a bound of ±10%.

$$Nu=0.556 R{e}^{0.90} P{r}^{0.333} R{R}^{0.26} {\epsilon }^{2.064} {\left(\frac{CC}{c}\right)}^{1.78} {\left(\frac{EC}{b}\right)}^{2.91} {\left(\frac{SC}{a}\right)}^{-0.315},R^2=0.9924$$

(29)

$$f=1.1687 R{e}^{-0.05656} P{r}^{-0.666} R{R}^{0.36} {\epsilon }^{0.50} {\left(\frac{CC}{c}\right)}^{0.752} {\left(\frac{EC}{b}\right)}^{0.4076} {\left(\frac{SC}{a}\right)}^{0.509}, R^2=0.9441$$

(30)

$$Nu=0.0334 R{e}^{0.87}{f}^{0.2215},R^2=0.9846$$

(31)

In general, Nu is strongly dependent on Re, with an exponent between 0.858 and 1.0 (Equations (23), (26) and (29)). Moreover, the friction factor decreases as Re increases with an exponent between −0.0673 and −0.0545 (Equations (24), (27) and (30)). The brick’s relative roughness influences the friction factor more than Nu, which is indicated by the exponent of the RR, where it is 0.36 and 0.26 in the friction factor and Nusselt number equations, respectively. Equations (29) and (30) show that the introduced geometrical groups ($\epsilon , \frac{CC}{c}, \frac{EC}{b}, \frac{SC}{a})$ proved to be significant in combining the data for the two settings into one equation. In particular, the friction factor correlation enhances considerably as its coefficient of determination increases substantially more than those for LDS and HDS.

Also, Nu is well correlated with the friction factor and Re, as given in Equation (31), with a high coefficient of determination (R² = 0.9846). Figure 7c shows the predictions of Nu against $R{e}^{0.87} f^{0.2215}$ with the line of slope 0.0529. The data points are distributed around the line in Figure 7c.

4.2. Effects of Design Parameters and Operating Conditions

4.2.1. Effect of Reynolds Number

The effect of Re on both Nu and f is investigated for both settings at Pr = 0.7, RR = 0.23, CC/c = 0.5172, and EC/b = 0.75. For LDS, SC/a = 1.3125, ε = 0.653, and Re = 184–10,200 while for HDS SC/a = 0.3175, ε = 0.48, and Re = 125–7600. Figure 8 shows the influence of the Reynolds number on the Nusselt number and the friction factor for both LDS and HDS. It is worth mentioning that the values of the Nusselt number and the friction factor are calculated from the general correlations (Equations (29) and (30)) using the parameters stated before. It can be noted that as the Reynolds number is increased, the friction factor tends to decrease, while the Nusselt number tends to increase. Figure 8a revealed that the effect of the Reynolds number on the Nusselt number is higher for the LDS than that for the HDS, especially at higher Reynolds numbers, while a minor difference between the Nusselt number of the HDS and LDS is observed at lower Reynolds numbers. The deviation between the LDS and HDS is due to the effect of voidage and stack channel in each setting. Figure 8b shows that the friction factor is relatively larger in the LDS than that in the HDS, although the pressure drop is much larger in the HDS case than that in the LDS. It can also be observed from Figure 8b that a significant effect of Reynolds number on the friction factor occurred at low Reynolds numbers while the effect is reduced as the Reynolds number is increased.

4.2.2. Effect of Brick’s Relative Roughness

The effect of RR on both Nu and f is investigated for both settings at Re = 5000, Pr = 0.7, CC/c = 0.5172, and EC/b = 0.75. For LDS, SC/a =1.3125 and ε = 0.653, while for HDS, SC/a = 0.3175 and ε = 0.48. Figure 9 shows the effect of the brick’s relative roughness on the Nusselt number and friction factor for both LDS and HDS. It is clear that as the brick’s relative roughness is increased, the Nusselt number and friction factor tend to increase. Moreover, the Nusselt number and the friction factor for LDS are larger than that for HDS. Figure 9 revealed that the slope of the Nusselt number is almost the same in both LDS and HDS cases, while the slope of the friction factor is lower in the HDS than that in the LDS. This indicates that the brick’s relative roughness has a slighter effect on the friction factor in the HDS case than in the LDS case.

4.2.3. Effect of the Stack Channel Group

The effect of SC/a on both Nu and f is investigated for both settings at Re = 5000, Pr = 0.7, RR = 0.50, CC/c = 0.5172, EC/b = 0.75, and SC/a = 0.3–1.3. ε = 0.653 for LDS and ε = 0.48 for HDS. Figure 10 shows the effect of the stack channel group on the Nusselt number and friction factor for both LDS and HDS. It is obvious that as the stack channel group is increased, the Nusselt number decreases while the friction factor tends to increase. Moreover, the values of the Nusselt number and the friction factor in the LDS are significantly larger than those in the HDS.

4.2.4. Effect of Tunnel Voidage

The effect of voidage on both Nu and f is investigated for both settings at Re = 5000, Pr = 0.7, CC/c = 0.5172, EC/b = 0.75, and ε = 0.45–0.70. For LDS, SC/a = 1.3125 and HDS SC/a = 0.3175. Figure 11 illustrates the dependency of the Nusselt number and friction factor for both the LDS and HDS upon the tunnel voidage. It can be noted from Figure 11 that as tunnel voidage is increased, both the Nusselt number and the friction factor increase.

4.2.5. Effect of Prandtl Number

Due to the limited variation of Pr, its effect on both Nu and f is investigated for both settings for two values of 0.68 and 0.73 at CC/c = 0.5172 and EC/b = 0.75. For LDS, SC/a = 0.3175 and ε = 0.48, while for HDS, SC/a = 1.3125 and ε = 0.653. Figure 12 displays the effect of the Prandtl number and Reynolds number on the Nusselt number and friction factor for both the LDS and HDS. It is easy to see that there is almost no effect of the Prandtl number on the Nusselt number, while the Nusselt number is directly proportional to the Reynolds number, as shown in Figure 12a,b. A remarkable decrease in the friction factor can be observed as the Prandtl number is increased for both LDS and HDS cases, as shown in Figure 12c,d. It can also be shown in Figure 12c,d that the magnitude of the reduction of the friction factor is more in the HDS compared with the LDS case. Furthermore, it can be seen from Figure 12c,d that the friction factor is inversely proportional to the Reynolds number for both LDS and HDS cases.

4.3. Comparison with the Literature Correlations

Figure 13 compares the developed correlations for Nu and f with available correlations in the literature. The generally developed correlations in the present work (Equations (29) and (30)) are plotted using the characteristic of LDS and HDS, as shown in Figure 13. The developed Nu correlation equation compares very well with Abou-Ziyan [10] and Vogt and Beckmann [25] correlations in the turbulent flow region, as shown in Figure 13a. While the Abou-Ziyan correlation (Re = 6000–33,000) lies between the LDS and HDS correlations, Vogt and Beckmann’s correlation almost coincides with the HDS correlation. It has to be stated that Vogt and Beckmann correlation (Equation (4)) for the turbulent flow (Re = 6100–50,700) around bricks on the channels side was translated as $\mathsf{\Phi }$ (in the present work) equals 2.19445. Thus, Equation (4) becomes:

$$Nu=\mathsf{\Phi }·0.018·R{e}^{0.8}=0.0395{Re}^{0.8}$$

(32)

In addition, Refaey et al. [17] correlations [Equations (7) and (8)] for the cooling zone, using a U-shape or side guide vanes, are compared after considering the effect of guide vanes = 1. Although Equations (7) and (8) are valid for Re between 13,609 and 27,634, they were slightly extrapolated between 10,000 and 11,000 to provide a way to compare them with the present work. Both correlations lie between the LDS and HDS of the present correlation.

Figure 13b shows the friction factor correlation of the present work for LDS and HDS compared against Abou-Ziyan’s correlation [10], which lies between the present data. It is worth mentioning that Equation (3) predicts negative values of friction factor for Re greater than or equal to 2500. It was stated by Vogt and Beckmann [24] that negative friction factors may be predicted for G₂ less than 0.5 (in the present work, G₂ = 0.4028). Therefore, it was not useful to plot it in Figure 13b.

In conclusion, the present correlations compare well with the literature correlations and extend the availability of data to low Reynolds numbers with simple forms to calculate Nu and f for brick settings in tunnel kilns.

5. Conclusions

The heat transfer and pressure drop data for tunnel kilns loaded with lattice bricks of different setting densities are developed and reported in dimensionless forms. The data are produced using three-dimensional CFD models for the two tested brick settings. The models are validated against experimental data reported in the literature. The developed correlations are useful for the design and operation of tunnel kilns to avoid kiln oversize and excessive specific energy consumption of the produced bricks. The following concluding remarks can be drawn from the present study:

• The developed Nusselt number and friction factor correlations considered the relative roughness of the bricks and the stack channels for the first time. Moreover, the correlations are valid for a practical range of Reynolds numbers that were not covered before in a simple and practical form.
• The data are correlated first for each setting density as a function of the Reynolds number, Prandtl number, and brick’s relative roughness. Then, the collected data for both settings were correlated with a single equation after introducing four geometrical parameters of the settings, including the setting voidage fraction and the ratios of column channel spacing, extension channel spacing, and stack channel spacing to brick length, width, and thickness, respectively. The general correlations are shown in Equations (29)–(31).
• The correlations are valid for Reynolds number between 125 and 10,200, Prandtl number between 0.68 and 0.73, brick’s relative roughness between 0.23 and 0.93, voidage fraction between 0.48 and 0.653, and the geometrical parameters of the tested lattice brick settings.
• The developed correlations of Nusselt numbers and friction factors are compared well with the available correlations in the literature in the valid range of parameters.
• In addition, the influence of the parameters considered in the developed correlations on Nusselt numbers and friction factors is investigated. It confirms that Nusselt numbers increase and the friction factors decrease significantly with Reynolds numbers and slightly with Prandtl numbers. At a constant Reynolds number, both the Nusselt number and the friction factor increase as the brick’s relative roughness is increased. On the other hand, as the stack channel group is increased, the Nusselt numbers decrease while the friction factors tend to increase. The voidage fraction of the setting has a monotonic effect on both Nusselt numbers and friction factors. Nusselt numbers for high-density are higher than those for low-density settings as the voidage fraction varies.