Numerical Study of Heated Tube Arrays in the Laminar Free Convection Heat Transfer

: Laminar free convection heat transfer from a heated cylinder and tube arrays is studied numerically to obtain the local and average Nusselt numbers. To verify the numerical simulations, the Nusselt numbers for a single cylinder were compared to other authors for the Rayleigh numbers of 10 3 and 10 4 . Furthermore, the vertically arranged heated tube arrays 4 × 1 and 4 × 2 with the tube ratio spacing S V /D = 2 were considered, and obtained average Nusselt numbers were compared to the existing correlating equations. A good agreement of the average Nusselt numbers for the single cylinder and the bottom tube of the 4 × 1 tube array is proved. On the other hand, the bottom tubes of the 4 × 2 tube array affect each other, and the Nusselt numbers have a different course compared to the single cylinder. The temperature fields for the tube array 4 × 4 in basic, concave, and convex configurations are studied, and new correlating equations were determined. The simulations were done for the Rayleigh numbers in the range of 1.3 × 10 4 to 3.7 × 10 4 with a tube ratio spacing S/D of 2, 2.5, and 3. On the basis of the results, the average Nusselt numbers increase with the Rayleigh numbers and tube spacing increasing. The average Nusselt number and total heat flux density for the convex configuration increase compared to the base one; on the other hand, the average Nusselt number decreases for the concave one. The results are applicable to the tube heaters constructional design in order to heat the ambient air effectively.


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 arrangement 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 0 ≤ Ra ≤ 10 7 .
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 affected by a bottom one, while changing Ra in the range of 2 × 10 4 -2 × 10 5 and the tube spacing from 2 to 9 cylinder diameters. It was found that the Nusselt numbers of the top cylinder are strongly affected by the tube spacing. A decrease of the Nusselt numbers occurs at small spacing, an enhancement prevails at higher one. Regarding the temperature differences, their effect 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 3 and 10 4 , 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 2 to 10 4 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 different 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 × 10 3 -1.6 × 10 5 .
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 centerto-center distance from 2 to more than 50 tube diameters were performed for Ra in the range between 5 × 10 2 and 5 × 10 5 . 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 2 and 10 4 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 3 and 3 × 10 3 . The Mach-Zehnder interferometer was used to visualize the temperature fields. It was found out that the location of each tube affects 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 effect of varying the Prandtl number and the vertical and horizontal pitch-to-diameter ratios for Ra of 1.5 × 10 8 . The heat transfer rates of the upper cylinder were affected by thermal plumes from the lower cylinders, and lower cylinders were not affected by the upper cylinders. When the vertical pitch is very small, the upper cylinders are affected.
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 coefficients 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. [12]. The vertical centerto-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 × 10 3 -5 × 10 5 . It was found out that the tube spacing affects 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 2 to 10 5 . 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 arrangement 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 3 and 10 4 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. 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 4 to 3.7 × 10 4 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' arrangement 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.

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. The mass (continuity equation), momentum equations, and the energy equation for incompressible and compressible flows are the following: Momentum equation in x direction: , Momentum equation in y direction: , Energy equation: .
The local Nusselt numbers Nuφ of any cylinder are calculated as [8]: . .
The average Nusselt numbers Nuav of any cylinder are calculated as [8]: 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]:

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) at the right symmetry line B-C (symmetry): b) on the cylinder surfaces (wall): c) at the bottom boundary line A-B (velocity inlet): d) at the left boundary line A-D (pressure outlet): e) at the top boundary line C-D (pressure outlet): 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 arrangement 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 sufficient 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. 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 differences 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.

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 3 and Ra = 10 4 . 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 3 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 4 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.
In comparison with Saitoh et al. [5], the discrepancies of Nuφ values are in the range of 0.02%-6.70% for Ra = 10 3 , and 0.36%-3.00% for Ra = 10 4 , 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 3 and 0%-5.13% for Ra = 10 4 , 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.

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 3 and Ra = 10 4 . 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 = 10 4 , Nuφ increases for all tubes I-IV in comparison with Ra = 10 3 . A significant increase of Nuφ on the downstream area of the tube II (φ = 0°-70°) is shown in comparison with Ra = 10 3 (Figure 5b). Due to the higher temperature difference between the tube surface and ambient air at Ra = 10 4 , 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 3 and 10 4 and the angles of φ = 0°-180°, Nuφ increases for Ra = 10 4 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. The ratio of the average Nusselt number for the i th cylinder and single cylinder / for a vertical column of four horizontal cylinders at Ra = 10 4 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 / , 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].
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 3 and 3.7% for Ra = 10 4 , respectively.

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 centerto-center dimensionless vertical distance SV/D = 2 and horizontal distance SH/D = 2 (Figure 1c)  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 3 are achieved on the tube I in the range of 1.167-3.615. For Ra = 10 4 , 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 = 10 3 . When comparing the same positions of the tubes at Ra = 10 3 and 10 4 and the angles of φ = 0°-360°, Nuφ increase for Ra = 10 4 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.
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: 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 = 10 3 . For Ra = 10 4 , 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 3 and Ra = 10 4 . The values of Nuav increase with Ra, specifically from 62% to 108% at Ra = 10 4 against Ra = 10 3 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 affect each other, where the values of Nuav are higher from 2% to 35% on each tube compared with the tube array 4 × 1. The temperature fields around the single cylinder and created thermal chimney at Ra = 10 3 and 10 4 are shown in Figure 8a 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 = 10 4 , 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 = 10 3 .

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 4 and 3.7 × 10 4 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. The ratio of presented Nu for the in-line bundle of horizontal cylinders and the single cylinder for Ra = 10 4 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 different 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 4 to 3.7 × 10 4 (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. 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 Figure 15 and Figure 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 affects 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 difference 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.
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, 15, and 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 4 and 3.7 × 10 4 is obvious in Table 4 and Table 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.
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 × 10 4 . 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 × 10 4 . 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 4 to 3.7 × 10 4 as follows: 1. The base configuration: with the percent standard deviation of error ESD = 0.70%, and error E ranged from -0.94% to +1.37%, 2. The concave configuration: with the percent standard deviation of error ESD = 0.65%, and error E ranged from -0.90% to +1.37%, 3. The convex configuration: with the percent standard deviation of error ESD = 0.54%, and error E ranged from -0.84% to +1.14%. 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.

Conclusions
The paper has reported on the results of the numerical investigation of laminar free convection heat transfer from various tube configurations. The single cylinder and the tube arrays of 4 × 1 and 4 × 2 arrangements at Ra = 10 3 and 10 4 were used for computational model verification. The main part of the research focused on the tube arrays in the configuration of 4 × 4 (base, concave, convex) at Ra in the range of 1.3 × 10 4 to 3.7 × 10 4 . From the obtained local and average Nusselt numbers and temperature fields, the following conclusions can be drawn: 1. The local Nusselt numbers Nuφ for the single circular cylinder at Ra = 10 3 and 10 4 were confirmed. 2. The average Nusselt numbers Nuav for the tube array 4 × 1 lie out of the error range, according to the Corcione correlating Equation [8], by 4.9% for Ra = 10 3 and 3.7% for Ra = 10 4 , respectively. 3. The average Nusselt numbers Nuav for the tube array 4 × 2 lie out of the error range, according to the Corcione correlating Equation [9], by 1.2% for Ra = 10 3 and lie in the error range for Ra = 10 4 . 4. There is a good agreement between the Nuφ of the single cylinder and bottom tube of the tube array 4 × 1 at Ra = 10 3 and 10 4 .

5.
The left bottom tube for the tube array 4 × 2 is affected by the right tube column; therefore, the Nuφ values have different courses for Ra = 10 3 and 10 4 . 6. The Nuav values increase with the S/D ratio and Ra increasing for the tube array 4 × 4 at the base, concave, and convex configurations. 7. The highest values of Nuav for the tube array 4 × 4 are noticed at the convex configuration which has discrepancies up to 3.7% in comparison with the concave one depending on Ra in the range of 1.3 × 10 4 to 3.7 × 10 4 . The value of Nuav for the base configuration lies between the convex and concave ones for all Ra. 8. The concave configuration achieves a decrease of Nuarray and qarray compared with the base one.
The convex configuration achieves an increase of these representative values compared with the base one. 9. The convex configuration achieves an increase of Nuarray and qarray compared with the concave one for all ratios S/D and Ra. When increasing the ratio S/D and Ra, a decrease of Nuarray and qarray percentage discrepancies between the individual configurations is reached. 10