**1. Introduction**

Laminar free convection heat transfer from heated tubes to air has been used in various technical applications such as heat exchangers, electronic components, heat storage equipment, and waste heat recovery systems. The advantage of free convection heat transfer is a creation of thermal comfort without a draft and dust fragments swirling in the space. This energy-saving way of heat transfer is reliable and the operation without fans is cost-effective and noise reducing. On the other hand, heat transfer at free convection is lower; therefore, a sufficiently large heat transfer surface and its suitable shape and arrangemen<sup>t</sup> is needed to fulfil effective heat transfer. Owing to this, the research of a chimney effect over the heat transfer surfaces is justified and required.

Several authors have dealt with the numerical calculations of heat transfer from a heated cylinder or tubes array at free convection. Churchill and Chu [1] developed an empirical expression for the average Nusselt number *Nuav* of a horizontal cylinder at presented Rayleigh numbers *Ra* and Prandtl numbers *Pr*. The temperature fields and airflows around the horizontal cylinder at free convection were researched by Morgan [2]. Kuehn and Goldstein [3] studied laminar free convection heat transfer from the horizontal cylinder by solving the Navier–Stokes and energy equations using an elliptical numerical procedure for 10<sup>0</sup> ≤ *Ra* ≤ 107.

The effect of vertical spacing on free convection heat transfer for a pair of heated horizontal tubes arranged one above the other was studied by Sparrow and Niethammer [4]. They investigated how the heat transfer characteristics of a top cylinder are a ffected by a bottom one, while changing *Ra* in the range of 2 × 104–2 × 10<sup>5</sup> and the tube spacing from 2 to 9 cylinder diameters. It was found that the Nusselt numbers of the top cylinder are strongly a ffected by the tube spacing. A decrease of the Nusselt numbers occurs at small spacing, an enhancement prevails at higher one. Regarding the temperature di fferences, their e ffect on the Nusselt numbers is higher at small tube spacing, unlike the higher one. Saitoh et al. [5] presented bench-mark solutions with accuracy to at least three decimal places for *Ra* = 10<sup>3</sup> and 104, where *Nu*φ, *Nuav*, isotherms, streamlines, and vorticities were obtained. The laminar free convection around an array of two isothermal tubes arranged above each other for *Ra* in the range of 10<sup>2</sup> to 10<sup>4</sup> was studied numerically by Chouikh et al. [6]. Decreased *Nu* at close spacing and enhanced *Nu* at a large one occurred at the upper tube. Within the same tube spacing, heat transfer at the upper tube increases with *Ra*. Herraez and Belda [7] used the holographic interferometry method for the temperature field visualization around the tubes of di fferent diameters when changing the surface temperature. Furthermore, the functions of an exponential form were defined, and *Nu* numbers were calculated for the range of *Gr*·*Pr* = 1.2 × 103–1.6 × 105.

Corcione [8] numerically studied a steady laminar free convection from horizontal isothermal tubes set in a vertical array. In this research, a computer code based on the SIMPLE-C algorithm was developed and simulations for the arrays of 2–6 tubes with the center-to-center distance from 2 to more than 50 tube diameters were performed for *Ra* in the range between 5 × 10<sup>2</sup> and 5 × 105. Moreover, the heat transfer correlating equations for individual tubes in the array and for the whole tube array were proposed. Furthermore, a steady laminar free convection from a pair of vertical tube arrays for *Ra* in the range of 10<sup>2</sup> and 10<sup>4</sup> was numerical studied by Corcione [9]. The pairs of tube arrays consisted of 1–4 tubes with the center-to-center horizontal and vertical spacing from 1.4 to 24 and from 2 to 12 tube diameters, respectively.

Ashjaee and Yousefi [10] experimentally investigated the free convection heat transfer from the horizontal tubes arranged above each other and the tubes shifted in a horizontal direction. The tube spacing varied from 2 to 5 tube diameter in a vertical direction and from 0 to 2 tube diameter in a horizontal direction for the inclined array between *Ra* = 10<sup>3</sup> and 3 × 103. The Mach–Zehnder interferometer was used to visualize the temperature fields. It was found out that the location of each tube a ffects the flow and heat transfer on the other tubes in the array. Heo and Chung [11] numerically investigated the natural convection heat transfer of two staggered cylinders for laminar flows. They used the Ansys Fluent software to examine the e ffect of varying the Prandtl number and the vertical and horizontal pitch-to-diameter ratios for Ra of 1.5 × 108. The heat transfer rates of the upper cylinder were a ffected by thermal plumes from the lower cylinders, and lower cylinders were not a ffected by the upper cylinders. When the vertical pitch is very small, the upper cylinders are a ffected.

Cernecky et al. [12] visualized the temperature fields around the two heated tubes arranged above each other at the surface temperature of 323 K. The local heat transfer coe fficients were determined from the holographic interferogram images of the temperature fields and the results were compared with numerical simulations. New criterion equations for calculating *Nu* of electrically heated tubes at free convection heat transfer were presented by Malcho et al. [13]. The vertical center-to-center distance of the tubes was gradually changed from 20 to 100 mm with a step of 20 mm. The surface temperatures of the tubes were 30, 45, 60, 75, 90, and 105 ◦C at a fluid temperature of 20 ◦C. Lu et al. [14] numerically investigated heat transfer from the vertical tube arrays (2–10 horizontal tubes) in molten salts at *Ra* = 2 × 103–5 × 105. It was found out that the tube spacing a ffects the average heat transfer rate around the whole tube array. The research resulted in a determination of the heat transfer dimensionless correlating equations for any individual tube in two vertically aligned horizontal tubes.

Kitamura et al. [15] experimentally investigated the free convection heat transfer around a vertical row of 10 heated tubes with the diameters of 8.4, 14.4, and 20.4 mm at modified *Ra* in the range of 5 × 10<sup>2</sup> to 105. The thermal plumes arising from the upstream tubes remained laminar throughout the rows when the gaps between the tubes were smaller than 20.6 mm. When the gaps were higher than 30.6 mm, the thermal plumes began to sway and undergo the turbulent transition on the halfway of the rows. Cernecky et al. [16] visualized the temperature fields in a set of four tubes arranged above each other at horizontal shift of 1/4 and 1/2 of the tube diameter by the holographic interferometry. Subsequently, the local heat transfer coefficients around the tubes were evaluated. The spacing between the tubes in vertical direction has a significant effect on the heat transfer parameters compared with the horizontal spacing of the tube centers. The heat transfer parameters on the other tubes in the direction of free convection flow varied depending on the tube position and geometry of the array. Ma et al. [17] studied free convection for a single cylinder by the Computational Fluid Dynamics with laminar and different turbulent models. Moreover, a thermal chimney was studied, and the effect of horizontal spacing arrangemen<sup>t</sup> of tubes on temperature and velocity fields was clarified.

The mentioned results were compared with those obtained from our simulations and are presented in this paper. First, a numerical simulation for a single cylinder at *Ra* = 10<sup>3</sup> and 10<sup>4</sup> was performed, and the local Nusselt numbers *Nu*φ were compared with other authors. Subsequently, the results of *Nuav* for the tube arrays 4 × 1 and 4 × 2 are presented according to arrangemen<sup>t</sup> in Figure 1. The main part of the presented contribution are the results of *Nu*φ and *Nuav* for the tube array 4 × 4 at basic, concave, and convex configurations. The simulations were performed for *Ra* in the range of 1.3 × 10<sup>4</sup> to 3.7 × 10<sup>4</sup> and for vertical and horizontal spacing between the tubes *S*/*D* = 2, 2.5 and 3. New correlating equations for the Nusselt numbers and the temperature fields for the tube arrays were created. The goal of this paper is to compare the base configuration of the tube arrays 4 × 4 (Figure 2a) with the concave (Figure 2b) and convex one (Figure 2c) and to increase the *Nuarray* and *qarray*, respectively, by means of a suitable configuration. Furthermore, the goal is to determine new correlating equations for the 4 × 4 tube array of basic, concave, and convex configurations with subsequent evaluation of the tubes' arrangemen<sup>t</sup> and spacing effects on *Nuarray*. The results are applicable to the design of heating equipment where free convection is used such as the tube heaters with a hot water circuit and electric heaters which are suitable for bathroom and toilet heating.

**Figure 1.** Arrangement of the single cylinder and tube arrays 4 × 1 and 4 × 2: (**a**) single cylinder; (**b**) 4 × 1 tube array; (**c**) 4 × 2 tube array.

**Figure 2.** Arrangement of the tube arrays 4 × 4: (**a**) base configuration; (**b**) concave configuration; (**c**) convex configuration.

#### **2. Mathematical Formulation**

Free convection heat transfer is conducted between the heated cylinder surface and ambient environment (air) and between the heated tube array and ambient air, respectively. In the mathematical formulation, the uniform and constant temperature of the cylinder surface (tube array surfaces) *Ts* is considered, and the ambient air is taken with the temperature *Tf* without the interfering airflow. The air, flown at free convection, is considered as steady, two-dimensional (*x, y*), laminar, and incompressible. The properties of air are shown in Table 1.

**Table 1.** The properties of air.


The mass (continuity equation), momentum equations, and the energy equation for incompressible and compressible flows are the following:

$$\frac{\partial(\rho u)}{\partial \mathbf{x}} + \frac{\partial(\rho v)}{\partial y} = 0,\tag{1}$$

Momentum equation in x direction:

$$
\rho \left( u \frac{\partial u}{\partial \mathbf{x}} + v \frac{\partial u}{\partial \mathbf{y}} \right) = \mu \left( \frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2} \right) - \frac{\partial p}{\partial \mathbf{x}} \tag{2}
$$

Momentum equation in y direction:

$$
\rho \left( \mu \frac{\partial v}{\partial \mathbf{x}} + v \frac{\partial v}{\partial y} \right) = \mu \left( \frac{\partial^2 v}{\partial \mathbf{x}^2} + \frac{\partial^2 v}{\partial y^2} \right) - \frac{\partial p}{\partial y} + \rho g \beta \left( T\_s - T\_f \right), \tag{3}
$$

Energy equation:

$$
\alpha \left( u \frac{\partial T}{\partial \mathbf{x}} + v \frac{\partial T}{\partial y} \right) = \alpha \left( \frac{\partial^2 T}{\partial \mathbf{x}^2} + \frac{\partial^2 T}{\partial y^2} \right). \tag{4}
$$

The local Nusselt numbers *Nu*φ of any cylinder are calculated as [8]:

$$Nu\_{\phi} = \frac{qD}{\lambda\_f \left(T\_s - T\_f\right)} = -\frac{\partial T}{\partial r}\Big|\_{r=0.5}.\tag{5}$$

The average Nusselt numbers *Nuav* of any cylinder are calculated as [8]:

$$Nu\_{\rm av} = \frac{Q}{\lambda\_f \pi \left(T\_s - T\_f\right)} = -\frac{1}{\pi} \int\_0^{\pi} \left.\frac{\partial T}{\partial r}\right|\_{r=0.5} d\phi. \tag{6}$$

The average Nusselt number of the whole array *Nuarray* is obtained as an arithmetic average value of *Nuav* from the individual tubes in the array [8]:

$$Nu\_{array} = \frac{1}{N} \sum\_{i=1}^{N} Nu\_{av} \,\tag{7}$$

#### **3. Arrangement of Investigated Tube Arrays and Boundary Conditions**

The numerical simulations of temperature fields and calculation of the local and average Nusselt numbers were performed in the Ansys Fluent software. The following boundary conditions are applied: (a)atthelineB-C

 right symmetry (symmetry):

$$\frac{\partial V}{\partial X} = 0; \; \mathcal{U} = 0; \; \frac{\partial T}{\partial X} = 0,\tag{8}$$

(b) on the cylinder surfaces (wall):

$$\mathcal{U} = 0; \; V = 0; \; T = 1,\tag{9}$$

(c) at the bottom boundary line A-B (velocity inlet):

$$\frac{\partial V}{\partial Y} = 0; \; \mathcal{U} = 0; \; T = 0 \quad \text{if} \quad V \ge 0,\tag{10}$$

(d) at the left boundary line A-D (pressure outlet):

$$V = 0; \ \frac{\partial \mathcal{U}}{\partial X} = 0; \ \frac{\partial T}{\partial X} = 0 \quad \text{if} \quad \mathcal{U} \ge 0,\tag{11}$$

(e) at the top boundary line C-D (pressure outlet):

$$\frac{\partial V}{\partial \boldsymbol{Y}} = 0; \; \mathcal{U} = 0; \; \frac{\partial T}{\partial \boldsymbol{Y}} = 0 \quad \text{if} \quad V \ge 0. \tag{12}$$

Inlet boundary: velocity inlet and constant temperature of 295 K. Outlet boundary, left wall: pressure outlet = 0 gauge pressure and constant temperature of 295 K. Right wall: symmetry. Tube walls: constant temperature from 313 to 373 K (Figures 1 and 2).

A preliminary study of the mesh grid size was carried out where quadrilateral elements were used. As shown in Table 2, five different numbers of elements were created around the half tube where the average Nusselt numbers of the single cylinder arrangemen<sup>t</sup> were being evaluated. On the basis of this study, 50 elements around the half tube are sufficient from the accuracy point of view. Owing to having more information about the local Nusselt numbers, 90 elements around the half tube (2◦

angle increment) were used despite the computational time increase. Furthermore, the element length, perpendicular to the tube, was studied. Table 2 shows the numbers of elements on the length of 8 mm perpendicular to the tube with 90 elements around the half tube. It is clear that 15 is a su fficient number of elements. Considering the elements quality (aspect ratio, angles, etc.), 20 elements were chosen. The average element length size close to the tube is 0.3 mm.


**Table 2.** Mesh independence test for single cylinder.

The governing Equations (1)–(4) considering the boundary conditions (8)–(12) were solved by the finite volume method using the Ansys Fluent software. Both transient and steady analyses were compared with each other from an accuracy point of view. The comparison was performed on the single cylinder simulation. While the transient analysis was done with variable air properties, the steady analysis was done with both variable and constant air properties. Regarding the variable properties, thermal conductivity and viscosity were considered as polynomial functions of temperature while specific heat was considered as a piecewise-polynomial function of temperature.

On the basis of negligible di fferences shown in Figure 3, the steady pressure-based analyses with the variable air properties were performed for the other arrangements. The Boussinesq approximation was used as a computational model and the pressure–velocity coupling was handled by the SIMPLE-C algorithm. For a momentum and energy gradient, the Quick and Second order upwind schemes were used. The solution was considered to be fully converged when the residuals of continuity, x-velocity, y-velocity, and energy met the convergence criterion 10−6. Finally, *Nu*φ were exported and subsequently evaluated.

**Figure 3.** The comparison of the local Nusselt numbers *Nu*φ of the single cylinder for various computational models at *Ra* = 10<sup>3</sup> and *Ra* = 104.

#### **4. Results and Discussion**

#### *4.1. Results of the Single Cylinder*

Laminar free convection heat transfer from the single isothermal circular cylinder with a diameter of 20 mm was studied (Figure 1a). The simulations were performed for *Ra* = 10<sup>3</sup> and *Ra* = 104. For a comparison, numerical simulations for the steady analysis with constant and variable material properties were performed (Figure 3).

The results of Corcione [8] were compared with ours. The discrepancies of *Nu*φ values for *Ra* = 10<sup>3</sup> are in the range of 0.37–5.90% for the steady state at constant air properties, and 0.14% to 6.11% for the steady state at variable air properties. The discrepancies of *Nu*φ values for *Ra* = 10<sup>4</sup> are in the range of 0.23–3.14% for the steady state at constant air properties and 0.33–3.07% for the steady state at variable air properties.

The present numerical simulations quantitatively match the data of the numerical simulations by Corcione [8], Saitoh et al. [5], and Wang et al. [18] for the single cylinder shown in Figure 4. The present numerical simulations were performed at the steady state analysis and variable air properties, described in Section 3.

**Figure 4.** The comparison of the local Nusselt numbers *Nu*φ of the single cylinder with other authors at *Ra* = 10<sup>3</sup> and *Ra* = 104.

In comparison with Saitoh et al. [5], the discrepancies of *Nu*φ values are in the range of 0.02–6.70% for *Ra* = 103, and 0.36–3.00% for *Ra* = 104, respectively. In comparison with Wang et al. [18], the discrepancies of *Nu*φ values are in the range of 1.05–7.84% for *Ra* = 10<sup>3</sup> and 0–5.13% for *Ra* = 104, respectively.

The numerical simulations for the single cylinder were performed to create a suitable computational model which was verified with other authors. On the basis of this model, the other simulations for the tube arrays 4 × 1, 4 × 2, and 4 × 4 were performed.

#### *4.2. Results of the Tube Array 4* × *1a*

We have studied the laminar free convection heat transfer from the vertical array of four isothermal circular tubes with a diameter of 20 mm. The simulations were performed for the tubes at a center-to-center dimensionless vertical distance of *SV*/*D* = 2 (Figure 1b) at *Ra* = 10<sup>3</sup> and *Ra* = 104. The courses of the local Nusselt numbers *Nu*φ for the individual tubes in the array are shown in Figure 5a,b for the angles of φ = 0◦ to 180◦. Owing to the model symmetry, the courses of *Nu*φ are identical for the right side as well.

As it is shown in Figure 5, the highest values of *Nu*φ are achieved on the tube I circumference. For *Ra* = 104, *Nu*φ increases for all tubes I–IV in comparison with *Ra* = 103. A significant increase of *Nu*φ on the downstream area of the tube II (φ = 0◦–70◦) is shown in comparison with *Ra* = 10<sup>3</sup> (Figure 5b). Due to the higher temperature difference between the tube surface and ambient air at *Ra* = 104, heat transfer to the tube II is affected by the heat flux from the tube I, and therefore *Nu*φ increases. Simultaneously, the tubes III and IV are affected by the heat flux of previous tubes. When comparing the same positions of the tubes at *Ra* = 10<sup>3</sup> and 10<sup>4</sup> and the angles of φ = 0◦–180◦, *Nu*φ increases for *Ra* = 10<sup>4</sup> in the range of 35–61% for the tube I, 30–171% for the tube II, 29–158% for the tube III, and 29–121% for the tube IV.

**Figure 5.** The course of the local Nusselt numbers *Nu*φ for the individual tubes in the tube array 4 × 1: (**a**) *Ra* = 103; (**b**) *Ra* = 104.

The ratio of the average Nusselt number for the *i*th cylinder and single cylinder *Nuiav*/*Nu*0*av* for a vertical column of four horizontal cylinders at *Ra* = 10<sup>4</sup> is shown in Figure 6a. The presented values were compared with Razzaghpanah and Sarunac [19] on the basis of Table 3 and Figure 9 [19] and Corcione [8] on the basis of Equation (17) [8]. The average *Nu* for the first cylinder in a column is the same as for the single cylinder. The relative error of the other presented tubes for the ratio *Nuiav*/*Nu*0*av*, compared with Razzaghpanah and Sarunac [19], is up to 5.2%. Figure 6b compares the presented average Nusselt numbers of the *i*th cylinder with Equation (22) of the authors Razzaghpanah and Sarunac [19] and Equations (11)–(13) of the authors Kitamura et al. [15].

**Figure 6.** Variation of the (**a**) ratio *Nuiav*/*Nu*0*av* vs. cylinder number; (**b**) average Nusselt number of the *i*th cylinder.

Corcione [8] obtained numerical results for the average Nusselt number of the tube array *Nuarray* which may be correlated to the Rayleigh number *Ra*, the tube ratio spacing *SV*/*D*, and the number of the tubes in the array *N*:

$$\begin{aligned} \text{Nu}\_{\text{array}} &= \text{Ra}^{0.235} \Big[ 0.292 \ln \Big| \left( \frac{S\_V}{D} \right)^{0.4} \times N^{-0.2} \Big] + 0.447 \Big], \\ &2 \le N \le 6; \quad 5 \times 10^2 \le \text{Ra} \le 5 \times 10^5, \\ &S\_V/D \le 10 - \log(\text{Ra}), \end{aligned} \tag{13}$$

with the percent standard deviation of error *ESD* = 2.25% and error *E* in the range of -4.79% to +5.27%. The values *Nuarray* obtained by our simulations lie out of the error range according to Corcione [8] by 4.9% for *Ra* = 10<sup>3</sup> and 3.7% for *Ra* = 104, respectively.

#### *4.3. Results of the Tube Array 4* × *2*

We have studied the laminar free convection heat transfer from the two vertical arrays of four isothermal circular tubes with a diameter of 20 mm. The simulations were performed for the center-to-center dimensionless vertical distance *SV*/*D* = 2 and horizontal distance *SH*/*D* = 2 (Figure 1c) at *Ra* = 10<sup>3</sup> and *Ra* = 104.

The courses of the local Nusselt numbers *Nu*φ for the individual tubes in the array 4 × 2 are shown in Figure 7. The courses of *Nu*φ are valid for the left tube array; the right tube array has mirrored *Nu*φ courses. As shown in Figure 7a, the highest values of *Nu*φ for *Ra* = 10<sup>3</sup> are achieved on the tube I in the range of 1.167–3.615. For *Ra* = 104, the local Nusselt numbers on the tube I are in the range of 1.518–6.118, which is an increase of 30–69% in comparison with *Ra* = 103. When comparing the same positions of the tubes at *Ra* = 10<sup>3</sup> and 10<sup>4</sup> and the angles of φ = 0◦–360◦, *Nu*φ increase for *Ra* = 10<sup>4</sup> in the range of 30–70% for the tube I, 48–183% for the tube II, 55–245% for the tube III, and 48–265% for the tube IV.

**Figure 7.** The course of the local Nusselt numbers *Nu*φ for the individual tubes in tube array 4 × 2: (**a**) *Ra* = 103; (**b**) *Ra* = 104.

Corcione [9] obtained numerical results for the average Nusselt number of the double tube array *Nuarray* which may be correlated to the Rayleigh number *Ra*, the tube ratio spacing *SV*/*D* and *SH*/*D*, and the number of the tubes in each vertical tube array of the pair *N*:

$$\begin{array}{ll} Nu\_{array} = 0.43 Ra^{0.235} (S\_H/D)^{0.14} (S\_V/D)^{0.2} \times N^{-0.1}, \\ 2 \le N \le 4; \; 2.4 - 0.2 \log (Ra) \le S\_H/D < 5, \\\ 2 \le S\_V/D < 5; \; 10^2 \le Ra \le 10^4. \end{array} \tag{14}$$

with the percent standard deviation of error *ESD* = 2.12% and error *E* in the range of -5.32% to +5.67%. The values *Nuarray* obtained by our simulations lie out of the error range according to Corcione [9] by 1.2% for *Ra* = 103. For *Ra* = 104, our values lie in the error range according to Corcione [9].

The values of average Nusselt numbers *Nuav* for the single cylinder and tube arrays 4 × 1 and 4 × 2 are given in Table 3 for *Ra* = 10<sup>3</sup> and *Ra* = 104. The values of *Nuav* increase with *Ra*, specifically from 62% to 108% at *Ra* = 10<sup>4</sup> against *Ra* = 10<sup>3</sup> for individual tubes within the same arrangement. For the tube array 4 × 1, the values of *Nuav* decrease from the tube I to the tube IV. For the tube array 4 × 2, the columns a ffect each other, where the values of *Nuav* are higher from 2% to 35% on each tube compared with the tube array 4 × 1.


**Table 3.** Average Nusselt numbers *Nuav* for the single cylinder and tube arrays 4 × 1, 4 × 2 at *Ra* = 10<sup>3</sup> and *Ra* = 104.

The temperature fields around the single cylinder and created thermal chimney at *Ra* = 10<sup>3</sup> and 10<sup>4</sup> are shown in Figure 8a,b. The temperature fields of the tube arrays 4 × 1 and 4 × 2 together with created thermal chimneys at *Ra* = 10<sup>3</sup> and 10<sup>4</sup> are shown in Figure 8c–f.

**Figure 8.** *Cont.*

**Figure 8.** The temperature fields of the cylinder and tube arrays 4 × 1, 4 × 2: (**a**) single cylinder at *Ra* = 103; (**b**) single cylinder at *Ra* = 104; (**c**) tube array 4 × 1 at *Ra* = 103; (**d**) tube array 4 × 1 at *Ra* = 104; (**e**) tube array 4 × 2 at *Ra* = 103; (**f**) tube array 4 × 2 at *Ra* = 104.

As the thermal chimney develops upwards, its temperature decreases due to heat transfer to the cold ambient air. A density difference between the thermal chimney and ambient air induces a buoyancy force; near the heated cylinder, the buoyancy force is higher.

In Figure 8c,d, it can be seen a gradual increase of the thermal boundary layer thickness in a vertical direction of the tube array 4 × 1. The widening of the thermal boundary layer causes the diminishing of the temperature gradient on the surfaces and the reduction of the local and average Nusselt numbers *Nu*φ and *Nuav*, respectively. In Figure 8e–f, a mutual effect of two tube columns in the tube array 4 × 2 can be seen, where an in-draft heat occurs between the first and second columns and a common thermal chimney is achieved over the tube array. For *Ra* = 104, a thinner thermal boundary layer around the tube circumferences can be observed. Owing to this, the values of *Nu*φ and *Nuav* increase compared with *Ra* = 103.

The comparison of the local Nusselt numbers *Nu*φ of the single cylinder with the bottom tube of the tube array 4 × 1 and with the bottom left tube of the tube array 4 × 2 is shown in Figure 9 for *Ra* = 10<sup>3</sup> and 104. The courses of *Nu*φ for the tube array 4 × 2 are valid for the left tube column; the right tube column has mirrored *Nu*φ courses. As shown, the values *Nu*φ for the bottom tube of the 4 × 1 array are similar to the values *Nu*φ for the single cylinder at *Ra* = 10<sup>3</sup> and 104. When comparing the bottom tube of the 4 × 2 array and the single cylinder, there are different courses due to the right tube column effect.

#### *4.4. Results of the Tube Array 4* × *4 (Base, Concave, Convex)*

The laminar free convection heat transfer from the four vertical arrays of four isothermal circular tubes with the diameter of 20 mm has been studied. Moreover, the effect of the tube ratio spacing on the Nusselt numbers had been investigated for different *Ra* before. The tube spacing in the vertical direction *SV* and horizontal one *SH* increased in the same way; therefore, general tube ratio spacing *S*/*D* is defined. The simulations were performed for the *S*/*D* ratios of 2, 2.5, and 3 (Figure 2) for the tube array 4 × 4 at three different configurations (Base—Figure 2a, Concave—Figure 2b, Convex—Figure 2c). The courses of the *Nu*φ values in the tube array 4 × 4 for base, concave, and convex configurations at *Ra* = 1.3 × 10<sup>4</sup> and 3.7 × 10<sup>4</sup> with a *S*/*D* ratio of 2 are shown in Figure 10. The courses for the same parameters and the *S*/*D* ratio of 2.5 are shown in Figure 11 and for the *S*/*D* ratio of 3 in Figure 12. When increasing the Rayleigh number, the convection increases in strength. Furthermore, the *Nuav* and *Nuarray* values increase with the *S*/*D* ratio and *Ra* increasing.

**Figure 9.** The comparison of the local Nusselt numbers *Nu*φ of the single cylinder with the bottom tube of the tube arrays at *Ra* = 10<sup>3</sup> and *Ra* = 104.

(**c**)

**Figure 10.** *Cont.*

**Figure 10.** The local Nusselt numbers *Nu*φ for different configurations of the tube array 4×4 with the tube ratio spacing S/D of 2: (**a**) base configuration at *Ra* = 1.3 × 104; (**b**) base configuration at *Ra* = 3.7 × 104; (**c**) concave configuration at *Ra* = 1.3 × 104; (**d**) concave configuration at *Ra* = 3.7 × 104; (**e**) convex configuration at *Ra* = 1.3 × 104; (**f**) convex configuration at *Ra* = 3.7 × 104.

**Figure 11.** *Cont.*

**Figure 11.** The local Nusselt numbers *Nu*φ for different configurations of the tube array 4×4 with the tube ratio spacing S/D of 2.5: (**a**) base configuration at *Ra* = 1.3 × 104; (**b**) base configuration at *Ra* = 3.7 × 104; (**c**) concave configuration at *Ra* = 1.3 × 104; (**d**) concave configuration at *Ra* = 3.7 × 104; (**e**) convex configuration at *Ra* = 1.3 × 104; (**f**) convex configuration at *Ra* = 3.7 × 104.

**Figure 12.** The local Nusselt numbers *Nu*φ for different configurations of the tube array 4 × 4 with the tube ratio spacing S/D of 3: (**a**) base configuration at *Ra* = 1.3 × 104; (**b**) base configuration at *Ra* = 3.7 × 104; (**c**) concave configuration at *Ra* = 1.3 × 104; (**d**) concave configuration at *Ra* = 3.7 × 104; (**e**) convex configuration at *Ra* = 1.3 × 104; (**f**) convex configuration at *Ra* = 3.7 × 104.

The ratio of presented *Nu* for the in-line bundle of horizontal cylinders and the single cylinder for *Ra* = 10<sup>4</sup> was compared with Razzaghpanah and Sarunac [20] on the basis of Figure 10 [20]. For the ratio *SL*/*D* and *ST*/*D* = 2, the values *Nuarray* = 4.726 and *Nuo* = 4.710 were obtained. The mentioned ratios *Nuarray*/*Nuo* = 1.003 and *Nuo*/*Nuarray* = 0.997 are close to those obtained by Razzaghpanah and Sarunac [20] (1.073 and 0.900, respectively).

Subsequently, the average Nusselt numbers of the whole array *Nuarray* for three di fferent configurations (base, concave, convex) with the tube ratio spacing *S*/*D* of 2, 2.5, and 3 were compared for *Ra* in the range of 1.3 × 10<sup>4</sup> to 3.7 × 10<sup>4</sup> (Figure 13). On the basis of the results, the *Nuarray* value increases with the increase of *Ra* and tube spacing. The highest values *Nuarray* are achieved at the convex configuration. The discrepancies of *Nuarray* between the concave and convex configurations are up to 3.7% depending on *Ra*. Owing to negligible error ranges for the correlating equations, these are interchangeable with similar discrepancies. The *Nuarray* values for the base configuration lie between the values for the concave and convex ones for all considered *Ra*.

**Figure 13.** The average Nusselt numbers *Nuarray* for the base, concave, and convex configurations of the tube array 4 × 4 with the di fferent tube ratio spacing *S*/*D*.

The temperature fields for the tube arrays 4 × 4 with the tube ratio spacing *S*/*D* of 2 obtained by the numerical simulations are shown in Figure 14. The temperature fields for the tube array 4 × 4 with the tube ratio spacing *S*/*D* of 2.5 and 3 are shown in Figures 15 and 16, respectively. For the middle columns of the tubes (V to VIII–Figure 2), it is clear that the thermal plume from the lower tubes a ffects the local and average Nusselt numbers of the upper tubes. The upper tubes create a quasi-forced convection which tends to a heat transfer increase from the upper tubes. It causes a lower temperature di fference between the tube surface and the incoming air, which tends to a heat transfer decrease at the tubes downstream that is of a major importance at close spacing.

**Figure 14.** The temperature fields of the tube arrays 4 × 4 with the tube ratio spacing *S*/*D* = 2: (**a**) base configuration at *Ra* = 1.3 × 104; (**b**) base configuration at *Ra* = 3.7 × 104; (**c**) concave configuration at *Ra* = 1.3 × 104; (**d**) concave configuration at *Ra* = 3.7 × 104; (**e**) convex configuration at *Ra* = 1.3 × 104; (**f**) convex configuration at *Ra* = 3.7 × 104.

**Figure 15.** The temperature fields of the tube arrays 4 × 4 with the tube ratio spacing *S*/*D* = 2.5: (**a**) base configuration at *Ra* = 1.3 × 104; (**b**) base configuration at *Ra* = 3.7 × 104; (**c**) concave configuration at *Ra* = 1.3 × 104; (**d**) concave configuration at *Ra* = 3.7 × 104; (**e**) convex configuration at *Ra* = 1.3 × 104; (**f**) convex configuration at *Ra* = 3.7 × 104.

**Figure 16.** The temperature fields of the tube arrays 4 × 4 with the tube ratio spacing *S*/*D* = 3: (**a**) base configuration at *Ra* = 1.3 × 104; (**b**) base configuration at *Ra* = 3.7 × 104; (**c**) concave configuration at *Ra* = 1.3 × 104; (**d**) concave configuration at *Ra* = 3.7 × 104; (**e**) convex configuration at *Ra* = 1.3 × 104; (**f**) convex configuration at *Ra* = 3.7 × 104.

The outer column's air flow (I to IV–Figure 2) is affected by the inner tube columns and the heat is directed to the central part of the tube array with a subsequent thermal plume creation. It may be seen that the Rayleigh number increases the "chimney effect". As the temperature fields around the tubes do not affect each other significantly, the boundary layer thickness decreases and *Nuav* increases with the higher ratio *S*/*D*. The effect of the tube spacing on the temperature field shapes is obvious in Figures 14–16. On the basis of the results, we can state that the higher tube spacing, the more effective heat transfer.

The effect of the *S*/*D* value and the tube array configuration for *Ra* = 1.3 × 10<sup>4</sup> and 3.7 × 10<sup>4</sup> is obvious in Tables 4 and 5. The *Nuav* value for the tubes I–VIII, the *Nuarray* value for whole tube arrays, and the *qarray* value were evaluated for all configurations and *S*/*D* ratios. When increasing the *S*/*D* ratio, also *Nuav*, *Nuarray*, and *qarray* increase, where the tube array configuration has a crucial effect.


**Table 4.** Average Nusselt numbers for different ratios *S*/*D* and configurations of 4 × 4 tube arrays at *Ra* = 1.3 × 104.

**Table 5.** Average Nusselt numbers for different ratios *S*/*D* and configurations of 4 × 4 tube arrays at *Ra* = 3.7 × 104.


As it is shown in Table 4, *Nuarray* and *qarray* values decrease by 1.27% at the concave configuration compared with the base one for *S*/*D* = 2 and *Ra* = 1.3 × 104. On the contrary, these representative values increase by 2.32% at the convex configuration compared with the base one for the same parameters. Additionally, they increase by 3.64% at the convex configuration compared with the concave one. For *S*/*D* = 2.5, the same comparison is carried out, where *Nuarray* and *qarray* decrease by 0.93% at the concave configuration and increase by 1.19% at the convex configuration compared with the base one. When comparing the convex configuration with the concave one, *Nuarray* and *qarray* increase by 2.14%. For *S*/*D* = 3, *Nuarray* and *qarray* decrease by 0.42% at the concave configuration and increase by 0.37% at the convex configuration compared with the base one. When comparing the convex configuration with the concave one, *Nuarray* and *qarray* increase by 0.79%.

As shown in Table 5, *Nuarray* and *qarray* values decrease by 1.51% at the concave configuration compared with the base one for *S*/*D* = 2 and *Ra* = 3.7 × 104. These representative values increase by 1.24% at the convex configuration compared with the base one. Additionally, they increase by 2.79% at the convex configuration compared with the concave one. For *S*/*D* = 2.5, *Nuarray* and *qarray* decrease by 0.77% at the concave configuration and increase by 0.63% at the convex configuration compared with the base one. When comparing the convex configuration with the concave one, *Nuarray* and *qarray* increase by 1.41%. For *S*/*D* = 3, *Nuarray* and *qarray* decrease by 0.36% at the concave configuration and increase by 0.39% at the convex configuration compared with the base one. When comparing the convex configuration with the concave one, *Nuarray* and *qarray* increase by 0.76%. It is obvious that the increase of the *S*/*D* ratio and *Ra* values causes a decrease of *Nuarray* and *qarray* percentage discrepancies between the individual configurations.

New correlating equations for *Nuarray* were created for all *Ra* values in the range of 1.3 × 10<sup>4</sup> to 3.7 × 10<sup>4</sup> as follows:

1. The base configuration:

$$Nu\_{array} = 0.475 \text{ } Ra^{0.2351} \text{ } (S/D)^{0.1907}\text{ },\tag{15}$$

with the percent standard deviation of error ESD = 0.70%, and error E ranged from -0.94% to +1.37%,

2. The concave configuration:

$$Nu\_{array} = 0.463 \, Ra^{0.2348} (S/D)^{0.2132} \, , \tag{16}$$

with the percent standard deviation of error ESD = 0.65%, and error E ranged from -0.90% to +1.37%,

3. The convex configuration:

$$Nu\_{array} = 0.520 \text{ Ra}^{0.2303} (S/D)^{0.1555},\tag{17}$$

with the percent standard deviation of error ESD = 0.54%, and error E ranged from -0.84% to +1.14%.

The correctness of the correlating equations is shown in Figure 17, where the values of *Nuarray* from the correlating equations and numerical simulations were compared.

Considering the correlating equations for the individual tubes in the array, there is an effort to create the equation containing the tube position in numerical order from bottom to top. For the Nusselt numbers of individual tubes situated in the outer columns, it is possible to assume an approximation by a linear function with the maximum discrepancy of 3% for specific tube position, spacing, configuration, and *Ra*. However, the monotony is not uniform throughout the tube's interval. The correlating equation indicates a decreasing of *Nuav* with increasing the tube position. On the contrary, in general, *Nuav* for the tube 3 is higher than the value for the tube 2. This statement is valid for all configurations, tube spacing and *Ra*. The Nusselt numbers of the tubes for the inner columns do not have the same course character when changing the parameters. On the basis of the previous facts, the correlating equations for *Nuav* are not created.

**Figure 17.** Comparison between presented correlating equations and numerical results.
