Natural convection in a cold enclosure with four hot inner cylinders based on diamond arrays (Part-I: Effect of horizontal and vertical equal distance of inner cylinders)

International Journal of Heat and Mass Transfer 111 (2017) 755–770

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Natural convection in a cold enclosure with four hot inner cylinders based on diamond arrays (Part-I: Effect of horizontal and vertical equal distance of inner cylinders) Gi Su Mun a, Yong Gap Park b, Hyun Sik Yoon c, Minsung Kim a, Man Yeong Ha a,⇑ a b c

School of Mechanical Engineering, Pusan National University, Jang Jeon 2-Dong, Geum Jeong Gu, Busan 609-735, Republic of Korea Rolls-Royce and Pusan National University Technology Centre in Thermal Management, Jang Jeon 2-Dong, Geum Jeong Gu, Busan 609-735, Republic of Korea Global Core Research Center for Ships and Offshore Plants, Pusan National University, Jang Jeon 2-Dong, Geum Jeong Gu, Busan 609-735, Republic of Korea

a r t i c l e

i n f o

Article history: Received 5 December 2016 Received in revised form 3 April 2017 Accepted 3 April 2017

Keywords: Natural convection Immersed boundary method Diamond array of cylinders

a b s t r a c t Two-dimensional numerical simulations were conducted to investigate the natural convection heat transfer induced by a temperature difference between the hot surfaces of four inner cylinders and the walls of an enclosure. The simulations were done using different Rayleigh numbers (103
Ra
106) and distances between neighboring cylinders (0.3
e
0.7). The immersed boundary method (IBM) was used to capture the virtual boundary of the four cylinders, which we arranged in a diamond array. Detailed analysis results including the distribution of streamlines, isotherms, and Nusselt numbers are presented. The thermal and flow fields reached steady state in the Rayleigh number range of 103
Ra
105, regardless of the variation in e. However, at Ra = 106, the numerical solutions show time-dependent characteristics when the distance between neighboring cylinders is 0.4
e
0.6. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Buoyancy-induced flow in an enclosure has a variety of applications, including heat exchangers, solar collector receivers, electronic equipment, and nuclear reactors. Many researchers have studied the effect of a single inner body inside the enclosure on the natural convection phenomena. The effects of various parameters on the natural convection have been instigated, such as the shapes and size of the inner body, the configuration of the enclosure, and the position of the inner body [1–13]. Many investigations have also examined the influence of multiple bodies in different positions on the natural convection characteristics within an enclosure. Lacroix and Joyeux [14] numerically studied the natural convection heat transfer for air from two vertically aligned, horizontal, heated cylinders inside a rectangular enclosure, which had vertical walls that had finite conductance and horizontal walls at the heat sink temperature in the Rayleigh number range of 103
Ra
106. They focused on the interaction between convection in the fluid-filled cavity and conduction in the vertical walls. Pelletier et al. [15], Persoons et al. [16], Reymond ⇑ Corresponding author. E-mail addresses: [email protected] (G.S. Mun), [email protected] (Y.G. Park), [email protected] (H.S. Yoon), [email protected] (M. Kim), [email protected] (M. Yeong Ha). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.04.004 0017-9310/Ó 2017 Elsevier Ltd. All rights reserved.

et al. [17], and Yoon et al. [18] investigated the natural convection phenomena in an enclosure with a pair of vertically aligned, horizontal, heated cylinders in different Rayleigh number ranges and considered the separation distance between the centers of two horizontal cylinders and their diameter. Park et al. [11] studied the characteristics of the natural convection with a temperature difference between two heated inner cylinders and a cold outer enclosure to reveal the effect of vertical changes in the position of two horizontally aligned, heated cylinders in the Rayleigh number range of 103
Ra
106. The bifurcation of natural convection from the steady state to the unsteady state depended on the Rayleigh number and the positions of the cylinders. Karimi et al. [22] numerically studied the effect on natural convection of the interaction from two heated horizontal cylinders to a cold square enclosure with focus on the dimensionless horizontal distance of the cylinders. They observed that the variation of the surface-averaged Nusselt number on the cylinder surface strongly depends on the distance between the cylinders. Other researchers focused on the effect of vertically aligned heated cylinders. Sadeghipour and Asheghi [24], Corcione [25], Kitamura et al. [26], and Ashjaee and Yousefi [27] investigated the natural convection heat transfer and fluid flow induced around vertical and inclined arrays of vertically aligned heated cylinders in different Rayleigh number ranges. The surface-averaged Nusselt number on the surface of the cylinders was strongly influenced

756

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

Nomenclature fi g L n Nu Nu hNui hNui P P Pr r R

momentum forcing acceleration of gravity [m=s2 ] length of square enclosure [m] normal direction to the wall local Nusselt number surface-averaged Nusselt number time-averaged local Nusselt number time- and surface-averaged Nusselt number pressure [Pa]  2 dimensionless pressure (¼ PqaL2 ) Prandtl number (¼ m=a) dimensionless radius of the cylinder (¼ R=L) radius of circular cylinder [m]

Ra t t T Th Tc ui ui xi

Rayleigh number (¼ gbL ðTmah T c Þ) time [s]  dimensionless time (¼ tL2a) dimensional temperature [K] hot temperature [K] cold temperature [K] velocity [m=s] u L dimensionless velocity (¼ ai ) Cartesian coordinates [m]

dimensionless Cartesian coordinates (¼ Li )

Greek symbols a thermal diffusivity [m2 =s] b thermal expansion coefficient [K 1 ] e dimensionless horizontal and vertical distance between the centers of four cylinders di2 Kronecker delta q density [kg=m3 ] m kinematic viscosity [m2 =s] u angle from the top of the circular cylinder c h dimensionless temperature (¼ TTT ) h T c

3

Subscripts/superscripts ⁄ dimensional value – surface-averaged quantity lower cyl lower cylinder middle right cyl middle right cylinder middle left cyl middle left cylinder upper cyl upper cylinder en enclosure

by the separation distance in the array. Park et al. [28] numerically studied the natural convection heat transfer induced by a temperature difference between a cold outer enclosure and four heated inner cylinders in a rectangular array with a Rayleigh number range of 103
Ra
106. They focused on the effects of the cylinder positions on the fluid flow and heat transfer in the enclosure. The flow and thermal fields became unstable depending on the cylinder position at high Rayleigh numbers of 105
Ra
106, in contrast to the cases at low Rayleigh numbers of 103
Ra
104. Few studies have investigated the effect of a diamond array of cylinders with different positions on the natural convection characteristics in an enclosure. Therefore, the present study investigates such a setup as continuation of research by Park et al. [28]. The focus is the results of the interaction between the multiple heated inner bodies. The fluid flow and heat transfer characteristics were analyzed by observing the isotherms and streamlines in the enclosure. The heat transfer capacity in terms of the Nusselt number was quantitatively estimated on the surfaces of the cylinders and enclosure walls.

2. Computation details 2.1. Numerical methods The immersed boundary method was applied to capture the virtual boundary of the cylinders in a Cartesian coordinate system, which is simpler and more efficient for dealing with a complex configuration compared to classical approaches such as bodyfitted curvilinear coordinate systems. We assume that the fluid is incompressible and Newtonian. The effects of thermal radiation, heat generation, and chemical reactions are neglected. Based on these assumptions, the non-dimensional governing equations of mass, momentum, and energy for an unsteady incompressible viscous flow and thermal fields are expressed as follows:

@ui q¼0 @xi

x

xi

ð1Þ

@ui @ui @P @ 2 ui ¼ þ Pr þ RaPrhðsin udi1 þ cos udi2 Þ þ f i þ uj @xi @t @xj @xj @xj ð2Þ @h @h @2h þ uj ¼ þh @t @xj @xj @xj

ð3Þ

The dimensionless variables are defined as follows:

t¼

t a L

2

; xi ¼

xi u L P L2 T  Tc ; h¼ ; ui ¼ i ; P ¼ L a qa2 Th  Tc

ð4Þ

In these equations,q, T, and a represent the density, dimensional temperature, and thermal diffusivity, respectively. The superscript ⁄ represents the dimensional variables. xi is the dimensionless Cartesian coordinate, ui is the corresponding dimensionless velocity component, t is the dimensionless time, P is the dimensionless pressure, and h is the dimensionless temperature. The nondimensionalization results in two dimensionless parameters, are Pr ¼ m=a and Ra ¼ gbL3 ðT h  T c Þ=ma, where m, g, and b are the kinematic viscosity, gravitational acceleration, and volume expansion coefficient, respectively. The terms q, f i , and h in Eqs. (1)–(3) are related to the immersed boundary method. The mass source/sink q in Eq. (1) and momentum forcing f i in Eq. (2) are imposed on the body surface and inside the body to satisfy the no-slip condition and mass conservation in the cell containing the virtual boundary. In Eq. (3), the heat source/sink h is applied to satisfy the isothermal boundary condition at the virtual boundary. The discretized equations and solution procedures are as follows:

^ ni  uni u 3 1 1 ^ ni Þ þ DIFðuni Þ Þ þ Pr½DIFðu ¼  NLðuni Þ þ NLðun1 i 2 2 2 Dt nþ1

þ RaPrhdi2 þ f i ^n @ 2 Pnþ1 @ u ¼ i  qnþ1 @xi @xi @xi

ð5Þ ð6Þ

757

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

^ ni  Dt unþ1 ¼u i

@Pnþ1 @xi

! ð7Þ

^hn  hn 3 1 1 nþ1 i i ¼  NLðhn Þ þ NLðhn1 Þ þ ½DIFðhnþ1 Þ þ DIFðhn Þ þ hi 2 2 2 Dt ð8Þ where uni is the intermediate time-step velocity, which is corrected , and superwith P nþ1 to satisfy mass conservation and becomes unþ1 i scripts n and n þ 1 represent the time levels. The nonlinear term NLðuÞ and diffusion term DIFðuÞ are defined as:

NLðui Þ ¼

@uj ui @xj

DIFðui Þ ¼

@ 2 ui @xj @xj

ð9Þ

A second-order linear or bilinear interpolation scheme was applied to satisfy the no-slip and isothermal conditions on the immersed boundary. Further details for the immersed boundary method are described by Kim et al. [31], Kim and Choi [32], and Choi [33]. A second-order accurate finite volume method was used for the spatial discretization of governing Eqs. (1)–(3). The fractional step method proposed by Choi and Moin [34] was used to simulate the time advancement of the flow field. In the discretization process, the advection terms were treated explicitly using the secondorder Adams-Bashforth scheme, and the diffusion terms were treated implicitly using the second-order accurate Crank-Nicolson scheme. Once the velocity and temperature fields are obtained, the local and surface-averaged Nusselt numbers are calculated using the equations below:

Nu ¼

 @h  @n

; Nu ¼ wall

1 S

Z

S

Nu ds

ð10Þ

0

where n represents the direction normal to the wall and S is the length of surface. 2.2. Computational conditions The computational domain, coordinate system, and the boundary conditions are shown in Fig. 1. The system consists of a square enclosure with four inner circular cylinders. The length of each wall of the enclosure is L, the radius of the cylinders is R = 0.1L, the Prandtl number is 0.7 (corresponding to that of air) and the

Fig. 2. A typical grid distribution at e = 0.3.

Rayleigh number range is 103
Ra
106. The positions of the four inner cylinders are changed along the horizontal and vertical directions based on the centerline of the enclosure in the range of  0.3
e
0.7, where e (¼ eL ) represents the dimensionless horizontal and vertical distances between the centers of four cylinders in each direction. The non-dimensionless temperatures are imposed on the cold walls (hc = 0) of the enclosure and surfaces (hh = 1) of the hot cylinders. No-slip and impermeability conditions were imposed on the surface of the cylinders and enclosure walls. All the fluid properties are assumed to be constant except for the density in the buoyancy term. The Boussinesq approximation was used to model the variation in the fluid density in the buoyancy term due to the change in the fluid temperature. Gravitational acceleration is applied in the negative y-direction. Fig. 2 shows the computational geometry in the x  y plane with a uniform grid distribution. A grid resolution of 501  501 along the horizontal (x) and vertical (y) directions was used. A grid dependency test was conducted with additional simulations on much finer grids with up to 601  601 points. Differences of less than 1% in the results of the surface-averaged Nusselt number were obtained using coarse and fine grids. 2.3. Validation test To validate the numerical method, preliminary numerical simulations were performed for the natural convection problem considered by Moukalled and Acharya [4] and Kim et al. [8]. Table 1 shows the comparison results in terms of the surface averaged

Table 1 Comparison of surface-averaged Nusselt numbers from the present study with previous numerical and experimental results. Ra

Mean Nusselt number at hot wall

3

Fig. 1. Computational domain, coordinate system, and along with boundary conditions.

10 104 105 106

Present study

Moukalled and Acharya [4]

Kim et al. [8]

5.107 5.128 7.836 14.462

– 5.49 8.377 15.4

5.093 5.108 7.767 14.11

758

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

Nusselt number on the cylinder surface in the Rayleigh number range of 103
Ra
106. Excellent agreement was achieved between the results. 3. Results 3.1. Steady and unsteady state flow regimes

Fig. 3. Bifurcation map state for different Rayleigh numbers and e.

Fig. 3 shows the bifurcation map in the Rayleigh number range of 103
Ra
106 and the dimensionless distance range of 0.3
e
0.7, which indicate a flow transition from the steady state to an unsteady state. The numerical solutions show timeindependent characteristics in the Rayleigh number range of 103
Ra
105 regardless of the variation in e. However, when the Rayleigh number increases to Ra = 106, the fluid flow accelerates due to the increase in the buoyancy effect. The numerical solutions show time-dependent characteristics at e = 0.4, e = 0.5, and e = 0.6 due to the combined effects of both the increase in the convection velocity and the dimensionless distances between cylinder centers.

(a) ε = 0.3

(f) ε = 0.3

(b) ε = 0.4

(g) ε = 0.4

(c) ε = 0.5

(h) ε = 0.5

(d) ε = 0.6

(i) ε = 0.6

(e) ε = 0.7

(j) ε = 0.7

Fig. 4. Distribution of isotherms and streamlines for different values of e at Ra = 103 (left column) and Ra = 104 (right column) (contour values range from 0 to 1 with 10 levels).

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

The fluid flow and thermal fields show the unsteady symmetric pattern at e = 0.6 and unsteady asymmetric pattern at e = 0.4 and e = 0.5 based on the vertical centerline of the enclosure at Ra = 106, which were obtained from the time-averaged field data reaching a statistically stationary state. 3.2. Fluid flow and thermal fields Fig. 4 shows the distribution of the isotherms and streamlines for different values of e in Ra = 103 and Ra = 104. A pair of primary vortices with two different vortex cores is formed on the left and right sides of the enclosure at e = 0.3, as shown in Fig. 4(a). The isotherms are evenly distributed and appear to have a symmetric pattern around the vertical centerline of the enclosure because conduction heat transfer is dominant rather than convection heat transfer. The two vortex cores within a pair of primary vortices become smaller when increasing from e = 0.3 to e = 0.5, but they become larger when increasing to e = 0.7. At e = 0.3, the distribution of the isotherms formed a radial shape that stretched out toward the center of the enclosure with increasing e. The thermal gradients between the surface of the cylinders and walls increased

759

due to the decrease in the distance between the cylinders and walls. Figs. 4(f)–(j) show that as the Rayleigh number increases to Ra = 104, the overall fluid flow and thermal structures in the natural convection phenomena are almost the same as those at Ra = 103. However, the buoyancy effect on the fluid flow and corresponding heat transfer slightly increase compared to that at Ra = 103. As a result, the thermal gradients near the lower part of the surface of the lower cylinder increases due to the increase in the buoyancy effect. Fig. 5 shows the distribution of the isotherms and streamlines for different values of e at Ra = 105 and Ra = 106. Figs. 5(a)–(e) show that when the Rayleigh number increases to 105, the magnitude of the convection velocity increases significantly due to the increase in the buoyancy effect compared to Ra = 103 and Ra = 104. At e = 0.3, a pair of primary vortices with two different vortex cores occurred on the left and right sides of the enclosure. The upper inner vortex is larger than the lower inner vortex due to the increase in the buoyancy effect. A single ascending plume is clearly visible above the surface of the upper cylinder due to the increase in the convection velocity in an enclosure.

(a) ε = 0.3

(f) ε = 0.3

(b) ε = 0.4

(g) ε = 0.4

(c) ε = 0.5

(h) ε = 0.5

(d) ε = 0.6

(i) ε = 0.6

(e) ε = 0.7

(j) ε = 0.7

Fig. 5. Distribution of isotherms and streamlines for different values of e at Ra = 105 (left column) and Ra = 106 (right column) (contour values range from 0 to 1 with 10 levels).

760

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

When e increases to e = 0.4 and e = 0.5, the distance between the surface of the upper cylinder and top wall of the enclosure becomes smaller. As a result, the two vortices at e = 0.4 divide into four small-scale secondary vortices at e = 0.5 near the top surface of the upper cylinder. This occurs to satisfy the mass continuity in the empty space caused by the flow separation of the primary vortices from the top surface of the upper cylinder. The two different vortex cores in each primary vortex shrink with increasing the e. The isotherms formed in the upper part of the enclosure are evenly distributed with increasing e due to the decrease in the distance between the surface of the upper cylinder and top enclosure wall. The distribution of the isotherms formed in the lower part of the enclosure move toward the center space of the cylinder array.

At e = 0.6, the small-scale secondary vortices above the top surface of the upper cylinder disappear due to the decrease in the distance from the top wall. When e increases further to e = 0.7, the two inner vortices in each primary vortex become larger because the fluid flow is confined near the corner spaces of the enclosure. This results from the increase in the space caused by the change in the cylinder positions. Compared to that at e = 0.6, the distribution of the isotherms formed in the upper part of the enclosure is more distorted due to the inner vortices located in each upper corner space of the enclosure. Figs. 5(f)–(j) show the distribution of the isotherms and streamlines for different values of e at Ra = 106. As shown in Fig. 3, the flow and thermal fields at e = 0.4, 0.5 and 0.6 for Ra = 106 show

Fig. 6. Distribution of time-averaged local Nusselt numbers along the surface of the upper cylinder (hNuiupper cyl ), the surface of the lower cylinder (hNuilower cyl ), the surface of the middle left cylinder (hNuimiddle left cyl ), the surface of the middle right cylinder (hNuimiddle right cyl ), and the walls of the enclosure (hNuien ) for different values of e at Ra = 103.

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

time-dependent characteristics. Therefore, the fluid flow and thermal fields with unsteady characteristics were obtained in the time-averaged field data after the time-averaging process was sufficiently achieved, as shown in Fig. 5(g)–(i). As a result, the asymmetric flow structures at e = 0.4 and e = 0.5 are attributed to the bifurcation problem. When the Rayleigh number increases to Ra = 106, the effect of convection on the fluid flow and the corresponding heat transfer intensify significantly compared to that at Ra = 105. This leads to a complex distribution of isotherms and streamlines. At e = 0.3, as shown in Fig. 5(f), a single ascending plume above the surface of the upper cylinder intensifies due to the increase in the convection velocity of the fluid circulating in the enclosure,

761

which leads to an increase in the thermal gradients near the middle of the top wall and the surfaces of the cylinders. A pair of two primary vortices with a single vortex core is formed on the left and right sides of the enclosure. The flow instability intensifies due to the combined effects of increased e and buoyancy, which leads to time-dependent characteristics for 0.4
e
0.6 at Ra = 106. The stable distribution of fluid flow in the enclosure at e = 0.3 becomes unstable as e increases from 0.3 to 0.4, resulting in the formation of asymmetric thermal and flow fields. The flow instability decreases with increasing e, which changes the flow characteristics from unsteady to steady. Therefore, the asymmetric thermal and flow fields in unsteady state at e = 0.4 and e = 0.5 change to a symmetric pattern in unsteady state when e = 0.6, which eventu-

Fig. 7. Distribution of time-averaged local Nusselt numbers along the surface of the upper cylinder (hNuiupper cyl ), the surface of the lower cylinder (hNuilower cyl ), the surface of the middle left cylinder (hNuimiddle left cyl ), the surface of the middle right cylinder (hNuimiddle right cyl ), and the walls of the enclosure (hNuien ) for different values of e at Ra = 104.

762

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

ally becomes a symmetric pattern in steady state as e increases to e = 0.7. At e = 0.4, two different primary vortices are identified in the enclosure, as shown in Fig. 5(g). The primary vortex occupying the right part of the enclosure shows a bi-cellular structure with two vortex cores and a single secondary vortex, which are separately positioned above the surface of the upper cylinder. Another primary vortex with a single vortex core filled the left part of the enclosure. Due to the distortion of the flow structure, an ascending thermal plume is generated above the upper cylinder toward the

upper-left part of the enclosure, which leads to thermal gradients concentrating near the middle-left of the top wall. At e = 0.5, the overall flow structure is similar to that at e = 0.4 except for a secondary vortex above the surface of the upper cylinder and the number of the inner vortex cores within a primary vortex that formed in the left part of the enclosure. However, the overall distribution of the isotherms is similar to that at e = 0.4. When e increases to e = 0.6, the asymmetric flow and thermal fields become symmetric patterns because the flow motion in an

Fig. 8. Distribution of time-averaged local Nusselt numbers along the surface of the upper cylinder (hNuiupper cyl ), the surface of the lower cylinder (hNuilower cyl ), the surface of the middle left cylinder (hNuimiddle left cyl ), the surface of the middle right cylinder (hNuimiddle right cyl ), and the walls of the enclosure (hNuien ) for different values of e at Ra = 105.

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

enclosure becomes stable relatively due to the positions of four cylinders closer to each wall. A pair of two primary vortices with two inner vortex cores is formed on the left and right sides of the enclosure. The upper inner vortex in each primary vortex shrinks compared to that at e = 0.5 due to the change in the locations of the cylinders. A pair of small secondary vortices is generated above the surface of the upper cylinder due to flow separation of the primary vortex from the upper surface of the upper cylinder. The distribution of isotherms is concentrated in

763

the gaps between each wall and cylinder surface due to the change in cylinder locations. When e increased further to e = 0.7, the overall flow structure and corresponding thermal field changed significantly. A pair of primary vortices with three inner vortex cores occurs on the left and right sides of the enclosure. The upper inner vortex becomes larger as e increases from 0.6 to 0.7 because a large amount of circulating fluid is tied up to near the upper side parts of the enclosure. Two inner vortices are located near the lower part of the

Fig. 9. Distribution of time-averaged local Nusselt numbers along the surface of the upper cylinder (hNuiupper cyl ), the surface of the lower cylinder (hNuilower cyl ), the surface of the middle left cylinder (hNuimiddle left cyl ), the surface of the middle right cylinder (hNuimiddle right cyl ), and the walls of the enclosure (hNuien ) for different values of e at Ra = 106.

764

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

middle left cylinder and middle right cylinder. As a result of the flow structure, the isotherms in the upper corner spaces are distorted and have a distribution with a rotating shape, which leads to a concentrated thermal gradient near the walls of the enclosure. The isotherms in the gaps between each wall of the enclosure and each cylinder surface are denser than those at e = 0.6. 3.3. Time-averaged local Nusselt number Fig. 6 shows the distributions of local Nusselt numbers along the surfaces of the upper cylinder (hNuiupper cyl), lower cylinder (hNuilower cyl), middle-left cylinder (hNuimiddle left cyl), middle-right cylinder (hNuimiddle right cyl), and enclosure walls (hNuien) for different values of e at Ra = 103. At Ra = 103, the isotherms are evenly distributed near the corner spaces of the enclosure, as shown in Fig. 4 (a)–(e). This indicates that the heat is predominantly transferred by the conduction rather than the convection. Therefore, the main factor for the distributions of the local Nusselt numbers is the distance between cylinder surfaces and walls of the enclosure in the horizontal and vertical directions. As a result, the values of hNuiupper cyl, hNuilower cyl, hNuileft cyl, hNuiright cyl, and hNu > en increase with increasing e. The thermal gradient in the inner space of the cylinders is smaller than that of the outer spaces of the enclosure, including the corner spaces. This is due to the confined inner space being surrounded by the cylinders. As shown in Fig. 6(a)–(d), low local Nusselt numbers occur around the inner space on the cylinder surfaces, regardless of the variations in e (u = 180° at the upper cylinder, u = 0° at the lower cylinder, u = 90° at the middle left cylinder,

(a)

and u = 270° at the middle right cylinder). However, the local Nusselt numbers increase in the outer part of the cylinders due to the smaller distance between the cylinder surfaces and the walls, regardless of the variations in e (u = 0° at the upper cylinder, u = 180° at the lower cylinder, u = 270° at the middle left cylinder, and u = 90° at the middle right cylinder). The maximum and minimum values of the local Nusselt number on the cylinder surfaces increase with e due to the increase in the thermal gradient in the space between the cylinder surfaces and the walls of the enclosure in the horizontal and vertical directions. High peak values of hNuien occur on the middle of the walls, regardless of the variations in e, because conduction is the dominant heat transfer mode rather than convection, and the values increase with e. Fig. 7 shows the distributions of local Nusselt numbers along the cylinder surfaces and walls for different values of e at Ra = 104. As the Rayleigh number increases to Ra = 104, the buoyancy effect on the fluid flow and corresponding heat transfer increased slightly compared to that at Ra = 103. This leads to increases in the local Nusselt number along the surfaces of the cylinders and walls, regardless of the variations in e. However, the heat is still predominantly transported by conduction rather than convection. As a result, the distance between the cylinder surfaces and walls in the horizontal and vertical directions is still the main factor. Therefore, the overall distributions of hNuiupper cyl, hNuilower cyl, hNuileft cyl, hNuiright cyl, and hNuien according to e at Ra = 104 are similar to those at Ra = 103. The local Nusselt numbers along the surfaces of the four cylinders increase slightly, regardless of the variations in the e compared to those at Ra = 103. This occurs due to the increase in convection velocity, as shown in Fig. 7(a)–(d).

(b)

(c) t = A

(d) t = B

(e) t = C

(f) t = D

Fig. 10. (a) Time histories, (b) power spectrum of the surface-averaged Nusselt number of the upper left cylinder, (c-f) instantaneous isotherms and streamlines at selected time instances at e = 0.4 and Ra = 106.

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

hNuien has single high and low peaks on the middle of the top wall and bottom wall of the enclosure, respectively, as shown in Fig. 7 (e). Fig. 8 shows the distributions of local Nusselt numbers along the cylinder surfaces and walls for different values of e at Ra = 105. The magnitude of the convection velocity increased significantly compared to those at Ra = 103 and Ra = 104, which changes the isotherm distribution and streamlines in the enclosure. Therefore, the distributions of the local Nusselt number for Ra = 105 are changed significantly compared Ra = 103 and Ra = 104. As shown in the distribution of hNuiupper cyl in Fig. 8(a), at e = 0.3, two high peaks of hNuiupper cyl occur around u = 70° and u = 290°. As e increases to 0.4, the maximum value of hNuiupper cyl occurs around u = 0° due to the increase in the thermal gradient near the top surface of the upper cylinder, which is caused by the downwelling plume above the surface of the upper cylinder. At e = 0.5 and e = 0.6, the distributions of hNuiupper cyl show similar patterns. The maximum value of hNuiupper cyl occurs around u = 0° due to the decrease in the distance between the surface of the upper cylinder and the top wall of the enclosure. The values of hNuiupper cyl distributed around 90°
u
270° increase compared to those at e = 0.3 and e = 0.4 because the heat transfer is accelerated on the lower surface of the upper cylinder due to the increase in e. As e increases further to e = 0.7, a single high peak of hNuiupper cyl occurs around u = 0°. This results from the heat transfer mainly occurring through the upper surface of the upper cylinder due to the very small distance from the top wall. This leads to a decrease in the value of hNuiupper cyl farther away from the top surface of the upper cylinder around 70°
u
290°.

(a)

765

As shown in the distribution of hNuilower cyl in Fig. 8(b), the maximum value of hNuilower cyl occurs around u = 180° regardless of the variations in e. The thermal gradient near the lower surface of the lower cylinder is larger than that of the upper surface due to the confined space surrounded by the heated cylinders. Therefore, the local Nusselt number on the surface of the lower cylinder sharply decreased toward the top surface from the middle of the side surface, regardless of the variations in e (around 0°
u
90° and 270°
u
360°). At e = 0.4, e = 0.5, and e = 0.6, two moderately high peaks of hNuilower cyl occur around u = 100° and u = 260°, in contrast to the cases of e = 0.3 and e = 0.7. As shown in the distribution of hNuimiddle left cyl in Fig. 8(c), at e = 0.3, the values of hNuimiddle left cyl are distributed around 0°
u
40° and 140°
u
360°. This occurs because the ascending flow motion of the primary vortex located on the left side of the enclosure generates a thinner thermal boundary layer near the left surface of the middle left cylinder. At e = 0.4 and e = 0.5, the overall distributions of hNuimiddle left cyl are similar. Two high peaks of hNuimiddle left cyl occur around u = 180° and u = 270° due to fluid impinging into the inner space of the cylinders with increasing e. At e = 0.6 and e = 0.7, a single high peak of hNuimiddle left cyl occurs around u = 270° due to the decrease in the distance between the surface of the middle left cylinder and the left wall. As shown in the distribution of hNuimiddle right cyl in Fig. 8(d), the distributions of hNuimiddle right cyl show the mirror-symmetry to those of hNuimiddle left cyl around the vertical centerline of the middle left cylinder, regardless of the variations in e. This results from the symmetric patterns of the isotherms and streamlines. As shown in the distribution of hNuien in Fig. 8(e), at e = 0.3, the maximum value of hNuien occurs on wall1 due to the presence of a

(b)

(c) t = A

(d) t = B

(e) t = C

(f) t = D

Fig. 11. (a) Time histories, (b) power spectrum of the surface-averaged Nusselt number of the upper left cylinder, (c-f) instantaneous isotherms and streamlines at selected time instances at e = 0.5 and Ra = 106.

766

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

single ascending plume above the upper cylinder. Two moderate high peaks of hNuien occur on wall2 and wall4 due to the presence of two inner vortices inside a pair of primary vortices in the enclosure. Low values of hNuien are distributed on wall3. At e = 0.4, two high peaks of hNuien occur on wall1 due to the secondary vortices formed near the middle of wall1. Even though the values of hNuien increase on wall2, wall3, and wall4 compared to those at e = 0.3, the distributions of hNuien show a similar pattern to those at e = 0.3 due to the similar flow structures near those walls in general. As e increases to e = 0.5, a moderately high peak of hNuien occurs on wall1 along with two high peaks, in contrast to e = 0.4. This results from the two secondary vortices above the upper cylinder which are divided into four secondary vortices due to the decrease in the distance between the surface of the upper cylinder and the top wall. The values of hNuien on wall2, wall3, and wall4 increase with similar patterns to those at e = 0.4. When e increases to e = 0.6, the distributions of hNuien change significantly compared to those in the range of 0.3
e
0.5. A single high peak of hNuien occurs on the middle of wall1. The positions of two moderately high peaks of hNuien on wall2 and wall4 are different from those in the range of 0.3
e
0.5, and the maximum value of hNuien on wall3 also increased noticeably. As e increases further to e = 0.7, the overall patterns of hNuien are similar to those at e = 0.6. However, due to the very small distance between the surfaces of the cylinders and each wall, the maximum values of hNuien on each wall increased significantly compared to those at e = 0.6. Fig. 9 shows the time-averaged local Nusselt numbers distributions at Ra = 106. As the Rayleigh number increases further up to Ra = 106, the distribution of the isotherms and streamlines changed

(a)

compared to those at lower Rayleigh numbers due to the increase in the buoyancy effect. As shown in the distribution of hNuiupper cyl in Fig. 9(a), at e = 0.3, two high peaks of hNuiupper cyl occur around u = 70° and u = 290°. As e increases to 0.4 and 0.5, a maximum value of hNuiupper cyl occurs around u = 180° because the heat transfer is accelerated on the lower surface of the cylinder due to the increases in e. The minimum value of hNuiupper cyl occurs around u = 300° due to the decrease of the thermal gradient near the cylinder surface caused by a single ascending plume that is slanted to the middle-left of the top wall from the upper cylinder surface. At e = 0.6, the values of hNuiupper cyl show mirror-symmetric patterns based on the vertical centerline of the enclosure, in contrast to the cases at e = 0.4 and e = 0.5. The values of hNuiupper cyl near the top surface of the upper cylinder are significantly increased compared to that at e = 0.5 due to the decrease in the gap between the top wall and the surface of the upper cylinder. Therefore, two high peaks of hNuiupper cyl occur around u = 0° and u = 180°. When e increases further up to 0.7, a single high peak of hNuiupper cyl occurs around u = 0°, and low values of hNuiupper cyl are distributed in the range of 70°
u
290°. As shown in the distribution of hNuilower cyl in Fig. 9(b), the maximum value of hNuilower cyl occurs around u = 180°, regardless of the variations in e. The value of hNuilower cyl sharply decreased toward the top surface from the middle of the side surface of the lower cylinder, regardless of the variations in e (around 0°
u
90° and 270°
u
360°), because the space near the upper part of the cylinder surface is surrounded by the heated cylinders. At e = 0.4, e = 0.5, e = 0.6, and e = 0.7, a single high peak of hNuilower cyl and two moderately high peaks of hNuilower cyl occur

(b)

(c) t = A

(d) t = B

(e) t = C

(f) t = D

Fig. 12. (a) Time histories, (b) power spectrum of the surface-averaged Nusselt number of the upper left cylinder, (c-f) instantaneous isotherms and streamlines at selected time instances at e = 0.6 and Ra = 106.

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

around u = 180° and around u = 100° and u = 260°, respectively. As shown in the distribution of hNuimiddle left cyl in Fig. 9(c), when e = 0.3, the values of hNuimiddle left cyl are distributed around 0°
u
40° and 140°
u
360° because the ascending flow motion of the primary vortex on the left hand side generates an thinner thermal boundary layer near the left surface of the middle left cylinder. At e = 0.4 and e = 0.5, the overall distributions of hNuimiddle left cyl are similar. A maximum value of hNuimiddle left cyl occurs around u = 180°. At e = 0.6, two high peak values of hNuimid-

767

dle left cyl occur around u = 180° and u = 290° due to the effect of two inner vortices within the primary vortex located in the left side of the enclosure. When e increases to e = 0.7, a single high peak of hNuimiddle left cyl occurs around u = 270° due to the decrease in the distance between the surface of the middle left cylinder and the left wall. As shown in the distribution of hNuimiddle right cyl in Fig. 9(d), The distributions of hNuimiddle right cyl show symmetric patterns to those of hNuimiddle left cyl based on the vertical centerline of the

Fig. 13. Time and surface-averaged Nusselt numbers as function of e for four different Rayleigh numbers at each wall: (a) top wall (hNuiT ), (b) bottom wall (hNuiB ), (c) right wall (hNuiR ), (d) left wall (hNuiL ), and (e) enclosure (hNuien ).

768

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

middle left cylinder, regardless of the variations in e, except for the cases of e = 0.4 and e = 0.5, which show asymmetric distribution of the isotherms and streamlines around the vertical centerline. At e = 0.4 and e = 0.5, minimum values of hNuimiddle right cyl occur around u = 340° and u = 320°, respectively. The maximum value of hNuimiddle right cyl occurs around u = 180°. As shown in the distribution of hNuien in Fig. 9(e), at e = 0.3, the maximum value of hNuien occurs on wall1 due to a single ascending plume above the upper cylinder. A single high peak of hNuien occurs on wall2 and wall4 near wall1 due to the increase in the buoyancy effect. Low values of hNuien are distributed on wall3. At e = 0.4 and e = 0.5, the high peaks of hNuien occur on wall1, and its values are located close to wall4 due to the ascending plume that is slanted to the middleleft of the top wall from the upper cylinder surface. The high peak of hNuien located on wall4 is larger than that at wall2. When e increases to e = 0.6, two high peaks of hNuien occur on wall1 due to two secondary vortices above the upper cylinder. High peaks of hNuien occur on wall2 and wall4 close to wall1. A single high peak and two moderately high peaks appear on wall3, in contrast to those at 0.3
e
0.5. This results from the decrease in the distance between the surface of the lower cylinder and the bottom wall. When e increases further to e = 0.7, the distribution of hNuien changes significantly compared to that at e = 0.6. Two high peaks of hNuien change to a single high peak and two moderately high peaks on wall1. A high peak of hNuien on wall3 also significantly increased compared to that at e = 0.6. A single high peak of hNuien is changed to a single high peak and moderately high peaks on wall2 and wall4 with the change in the flow structure. The maximum values are located on the middle of each wall.

3.4. Unsteady characteristics at Ra=106 Figs. 10–12 show the time histories and power spectrum of the surface-averaged Nusselt number of the upper cylinder (Nuupper cyl ), as well as the instantaneous isotherms and streamlines at selected time instances for 0.4
e
0.6 in Ra = 106. At e = 0.4, the time histories of Nuupper cyl oscillate periodically as a function of time, as shown in Fig. 10(a). The power spectrum reveals a primary frequency of approximately 45 and two different harmonics with a low power density near the primary frequency, as shown in Fig. 10(b). The instantaneous flow and thermal fields at times A, B, C, and D in Fig. 10(a) are shown in Fig. 10(c)–(f). The values of Nuupper cyl change as a function of time due to the oscillating motion of plumes above the surface of the upper cylinder. At t = A, Nuupper cyl has a maximum value due to increases in the thermal gradient on the cylinder surface caused by the downwelling plume above the surface of the upper cylinder. At t = B, the inner vortices within the primary vortex on the right side of the enclosure are shifted towards the top right corner of the enclosure when compared to those at t = A, resulting in a decrease in the thermal gradient on the surface of the upper cylinder. Therefore, the value of Nuupper cyl significantly decreases when t increases from A to B. At t = C, an inner vortex is located above the upper cylinder, which leads to an increase in Nuupper cyl . The value decreases when increasing from t = C to t = D due to the decrease in the thermal gradient on the surface of the upper cylinder caused by the small-scale secondary vortices above the upper cylinder. At e = 0.5, the time histories of Nuupper cyl oscillate periodically as a function of time, which is similar to those at e = 0.4, as shown in

Fig. 14. Time and surface-averaged Nusselt numbers as function of e for four different Rayleigh numbers: (a) upper cylinder (hNuiupper cyl ), (b) lower cylinder (hNuilower cyl ), (c) left cylinder (hNuimiddle left cyl ), and (d) right cylinder (hNuimiddle right cyl ).

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

Fig. 11(a). The power spectrum reveals primary and secondary frequencies of approximately 73 and 32, respectively, with harmonics and a very low power density at frequencies close to the primary frequency, as shown in Fig. 11(b). The instantaneous flow and thermal fields at times A, B, C, and D in Fig. 11(a) are shown in Fig. 11 (c)–(f). The time histories of Nuupper cyl at e = 0.5 show similar patterns to those at e = 0.4. Therefore, the overall flow structures at t = A and B show similar distributions to those at e = 0.4, as shown in Fig. 11(c) and (d). From t = B to t = C, the inner vortex in the right part of the upper cylinder disappears, which leads to the formation of a vortex circulating near the middle-right cylinder with two inner vortices. Therefore, the heat transfer on the surface of the upper cylinder is accelerated, resulting in an increase of Nuupper cyl from t = B to t = C. When t increases further to D, the vortex circulating near the middle-right cylinder disappears due to the large inner vortex to the right of the upper cylinder, which leads to a decrease in Nuupper cyl . At e = 0.6, the power spectrum shows a primary frequency of approximately 133. Therefore, oscillating motion as a function of time is modulated from a sinusoidal wave pattern with a primary frequency, as shown in Fig. 12(a) and (b). The instantaneous flow and thermal fields at times A, B, C, and D in Fig. 12(a) are shown in Figs. 12(c)–(f). Nuupper cyl has a maximum value at t = A. When t increases from A to C, as shown in Fig. 12(c)–(e), the convection velocity moving along the surface of the upper cylinder decreases due to the change in the vortex structure near the cylinders. Therefore, the values of Nuupper cyl decrease from t = A to t = C. When t increases to D, the inner vortices generated near the left and right parts of the upper cylinder at t = C disappear, resulting in acceleration of the heat transfer on the upper cylinder surface. 3.5. Time- and surface-averaged Nusselt number Fig. 13 shows the time- and surface-averaged Nusselt number on the top wall (hNuiT ), bottom wall (hNuiB ), right wall (hNuiR ), left wall (hNuiL ), and the entire enclosure (hNuiEN ) as a function of e in the Rayleigh number range of 103
Ra
106. The surfaceaveraged Nusselt numbers generally increased with the Rayleigh number, regardless of the walls and variation in e except for hNuiB . This was due to the acceleration in heat transfer caused by the increased convection velocity. As shown in the distribution of hNuiT in Fig. 13(a), at Ra = 103 and Ra = 104, the values of hNuiT show similar distributions with the variation in e. However, the values of hNuiT increase slightly with increasing Rayleigh number from Ra = 103 to Ra = 104, regardless of the variation in e due to the increase in convection velocity in the enclosure. At Ra = 105 for 0.3
e
0.5, the values of hNuiT are almost the same, regardless of the variation in e due to the presence of the secondary vortices positioned above the upper cylinder, which leads to a decrease in the heat transfer through the top wall. When e increases to 0.6 and 0.7, the decrease in the distance between the surface of the upper cylinder and the top wall leads to the secondary vortices above the upper cylinder disappearing, which increases hNuiT with increasing e. At Ra = 106, the distribution of hNuiT is significantly different from those at lower Rayleigh number. At e = 0.3, hNuiT is high on the top wall due to the ascending plume above the surface of the upper cylinder. When e increases from 0.3 to 0.4, there is a strong ascending plume stretching to the middle-left of the top wall caused by bifurcation from the increased flow instability in the enclosure. Thus, hNuiT increased slightly. At e = 0.5 and e = 0.6, the plume above the upper cylinder becomes weak, resulting in a decrease of hNuiT . When e increases

769

from to 0.6 and 0.7, hNuiT increases significantly compared to those at e = 0.5 and e = 0.6. As shown in the distribution of hNuiB in Fig. 13 (b), the value of hNuiB increases with increasing e and generally decrease with increasing Rayleigh number due to the increase in the buoyancy effect. The overall distributions of hNuiB for different Rayleigh numbers are similar, regardless of the variation in e except for the range of 0.3
e
0.5 in the Rayleigh number range of 105
Ra
106. As shown in the distributions of hNuiR and hNuiL in Fig. 13 (c) and (d), the convection velocity of the fluid moving along the side walls increases, so hNuiR and hNuiL increase with increasing Rayleigh number. They also increased with increasing e except for the range of 0.4
e
0.5 at Ra = 106, which results from the thermal and flow fields slanted to the left side. As a result of the distributions of the surface- and time-averaged Nusselt number on each wall, hNuiEN has a similar trend for each wall of increasing with the Rayleigh number and e except for the range of 0.4
e
0.6 at Ra = 106, as shown in the distribution of hNuiEN in Fig. 13(e). Fig. 14 shows the time- and surface-averaged Nusselt number on the cylinder surfaces as a function of e for the Rayleigh number range of 103
Ra
106. The surface-averaged Nusselt numbers generally increased with increasing Rayleigh number, regardless of the position of the cylinders and the variation in e. This was due to the acceleration in heat transfer caused by the increase of convection velocity except for those of the upper cylinder. As shown in the distribution of hNuiupper cyl in Fig. 14(a), hNuiupper cyl decreases from Ra = 103 to Ra = 105 due to the decrease in thermal gradient near the surface of the upper cylinder, which is caused by the interaction of an ascending plume from the lower cylinder with the upper cylinder, in addition to the position of the cylinders. The presence of an ascending plume above the lower cylinder generally reduces the heat transfer of the upper cylinder compared to that of the lower cylinder. However, when Ra = 106, the high convection velocity circulating along the surface of the upper cylinder increases, so the values of hNuiupper cyl increase compared to those at 103
Ra
106. hNuiupper cyl increases with increasing e, regardless of the variation in the Rayleigh number. As shown in the distribution of hNuilower cyl in Fig. 14(b), hNuilower cyl increases with increasing Rayleigh number and e except for the range of 0.4
e
0.5 at Ra = 106. As shown in Fig. 14(c) and (d), the distributions of hNuimiddle left cyl and hNuimiddle right cyl are similar, regardless of the variation in the Rayleigh numbers and e except for the range of 0.4
e
0.5 at Ra = 106 due to the symmetric patterns of the thermal and flow fields around the vertical centerline of the enclosure. hNuimiddle left cyl and hNuimiddle right cyl generally increased with increasing Rayleigh number and e. However, at 0.4
e
0.6 and Ra = 106, hNuimiddle left cyl and hNuimiddle right cyl gradually decrease with increasing e due to the p inner vortices near the middle left and right cylinders.

4. Conclusions In the present study, a two-dimensional numerical analysis was carried out to examine the natural convection induced by hot circular cylinders inside a cold square enclosure arranged in a diamond array for the Rayleigh number range of 103
Ra
106. The immersed boundary method based on the finite volume method was used to capture the virtual boundary of the cylinders in the Cartesian coordinate system. The cylinders were located on the vertical and horizontal centerlines of the enclosure, and they were

770

G.S. Mun et al. / International Journal of Heat and Mass Transfer 111 (2017) 755–770

moved along the centerlines to vary the distance between them in the range of 0.3
e
0.7. The thermal and flow fields changed from steady state to unsteady state depending on the Rayleigh number and e. The thermal and flow fields show the steady state solutions and symmetric structures based on the vertical centerline of the enclosure, regardless of the variation in the Rayleigh number and e except for 0.4
e
0.6 at Ra = 106. At e = 0.3 and Ra = 106, the numerical solution shows the time-independent characteristics, which changed to an unsteady state at 0.4
e
0.6. This resulted from the increase in the flow instability caused by the combined effect of the distance between neighboring cylinders and increased buoyancy. The unsteady characteristics returned to the steady state at e = 0.7. The thermal and flow fields are dependent on the Rayleigh number and e. At Ra = 103 and Ra = 104, the distribution of the isotherms and streamlines mainly depend on the distance between the cylinders and the cold walls because the conduction is a dominant heat transfer mode rather than the convection at the relatively low Rayleigh numbers. However, at Ra = 105 and Ra = 106, the distribution of the isotherms and streamlines are strongly dependent on the combined effects of the distance between neighboring cylinders and increased buoyancy in addition to the mutual interaction between the cylinders, including the interaction between the thermal plumes above the lower and upper cylinders. As a result, the distributions of isotherms and streamlines at Ra = 105 and Ra = 106 are different from those at Ra = 103 and Ra = 104. This leads to the distribution of the local Nusselt number along the cylinder surfaces and the walls depending on the Rayleigh number and e. The time and surface-averaged Nusselt number on walls and cylinder surfaces generally increased with increasing e due to the increase in the effect of the distance between the cylinder surfaces and the walls. The numerical solutions show time-independent characteristics except for the cases of 0.4
e
0.6 at Ra = 106. As a result, the time- and surface-averaged Nusselt number on walls and cylinder surfaces varies with different patterns in this range compared to 103
Ra
105 due to the combined effects of the distance between neighboring cylinders and increased buoyancy. Acknowledgement This subject is supported by Korea Ministry of Environment (MOE) as ‘‘the Chemical Accident Prevention Technology Development Project” (No. 2015001950002).

[8] B.S. Kim, D.S. Lee, M.Y. Ha, H.S. Yoon, A numerical study of natural convection in a square enclosure with a circular cylinder at different vertical locations, Int. J. Heat Mass Transfer 51 (2008) 1888–1906. [9] D.H. Kang, M.Y. Ha, H.S. Yoon, C. Choi, Bifurcation to unsteady natural convection in square enclosure with a circular cylinder at Rayleigh number of d, Int. J. Heat Mass Transfer 64 (2013) 926–944. [10] J.M. Lee, M.Y. Ha, H.S. Yoon, Natural convection in a square enclosure with a circular cylinder at different horizontal and diagonal locations, Int. J. Heat Mass Transfer 53 (2010) 5905–5919. [11] Y.G. Park, M.Y. Ha, H.S. Yoon, Study on natural convection in a cold square enclosure with a pair of hot horizontal cylinders positioned at different vertical locations, Int. J. Heat Mass Transfer 65 (2013) 696–712. [12] H.S. Yoon, M.Y. Ha, B.S. Kim, D.H. Yu, Effect of the position of a circular cylinder in a square enclosure on natural convection at Rayleigh number of Ra ¼ 105 , Phys. Fluids 21 (2009) 047101. [13] H.J. Lee, J.H. Doo, M.Y. Ha, H.S. Yoon, Effects of thermal boundary conditions on natural convection in a square enclosure with an inner circular cylinder locally heated from the bottom wall, Int. J. Heat Mass Transfer 65 (2013) 435–450. [14] M. Lacroix, A. Joyeux, Coupling of wall conduction with natural convection from heated cylinders in a rectangular enclosure, Int. Comm. Heat Mass Transfer 23 (1996) 143–151. [15] Q. Pelletier, D.B. Murray, T. Persoons, Unsteady natural convection heat transfer from a pair of vertically aligned horizontal cylinders, Int. J. Heat Mass Transfer 95 (2016) 693–708. [16] T. Persoons, I.M. O’Gorman, D.B. Donoghue, G. Byrne, D.B. Murray, Natural convection heat transfer and fluid dynamics for a pair of vertically aligned isothermal horizontal cylinders, Int. J. Heat Mass Transfer 54 (2011) 5163– 5172. [17] O. Reymond, D.B. Murray, T.S. O’Donovan, Natural convection heat transfer from two horizontal cylinders, Exp. Therm. Fluid Sci. 32 (2008) 1702–1709. [18] H.S. Yoon, J.H. Jung, Y.G. Park, Natural convection in a square enclosure with two horizontal cylinders, Numerical Heat Transfer 62 (2012) 701–721. [22] F. Karimi, H.T. Xu, Z. Wang, M. Yang, Y. Zhang, Numerical simulation of unsteady natural convection from heated horizontal circular cylinders in a square enclosure, Numerical Heat Transfer 65 (2014) 715–731. [24] M. Sadeghipour, M. Asheghi, Free convection heat transfer from arrays of vertically separated horizontal cylinders at low Rayleigh numbers, Int. J. Heat Mass Transfer 37 (1994) 103–109. [25] M. Corcione, Correlating equations for free convection heat transfer from horizontal isothermal cylinders set in a vertical array, Int. J. Heat Mass Transfer 48 (2005) 3660–3673. [26] K. Kitamura, A. Mitsuishi, T. Suzuki, F. Kimura, Fluid flow and heat transfer of natural convection induced around a vertical row of heated horizontal cylinders, Int. J. Heat Mass Transfer 92 (2016) 414–429. [27] M. Ashjaee, T. Yousefi, Experimental study of free convection heat transfer from horizontal isothermal cylinders arranged in vertical and inclined arrays, Heat Transfer Eng. 28 (2007) 460–471. [28] Y.G. Park, M.Y. Ha, J. Park, Natural convection in a square enclosure with four circular cylinders positioned at different rectangular locations, Int. J. Heat Mass Transfer 81 (2015) 490–511. [31] J. Kim, D. Kim, H. Choi, An immersed-boundary finite volume method for simulations of flow in complex geometries, J. Comp. Phys. 171 (2001) 132– 150. [32] J. Kim, H. Choi, An immersed-boundary finite-volume method for simulation of heat transfer in complex geometries, KSME Int. J. 18 (2004) 1026–1035. [33] C. Choi, A Numerical Study of the Effect of a Circular Cylinder on Natural Convection in Rectangular Channel Heated from Bottom Wall, Ph.D thesis, Pusan National University, (2014) pp. 1–171. [34] H. Choi, P. Moin, Effects of the computational time step on numerical solutions of turbulent flow, J. Comput. Phys. 113 (1994) 1–4.

References Further reading [1] J.R. Lee, M.Y. Ha, A numerical study of natural convection in a horizontal enclosure with a conducting body, Int. J. Heat Mass Transfer 48 (2005) 3308– 3318. [2] C. Shu, Y.D. Zhu, Efficient computation of natural convection in a concentric annulus between an outer square cylinder and an inner circular cylinder, Int. J. Numer. Meth. Fluids 38 (2002) 429–445. [3] D. Angeli, P. Levoni, G. Barozzi, Numerical predictions for stable buoyant regimes within a square cavity containing a heated horizontal cylinder, Int. J. Heat Mass Transfer 51 (2008) 553–565. [4] F. Moukalled, S. Acharya, Natural convection in the annulus between concentric horizontal circular and square cylinders, J. Thermophys. Heat Transfer 10 (1996) 524–531. [5] X. Xu, Z. Yu, Y. Hu, L. Fan, K. Cen, A numerical study of laminar natural convective heat transfer around a horizontal cylinder inside a concentric airfilled triangular enclosure, Int. J. Heat Mass Transfer 53 (2010) 345–355. [6] A. Kumar De, A. Dalal, A numerical study of natural convection around a square, horizontal, heated cylinder placed in an enclosure, Int. J. Heat Mass Transfer, Available online (2006). [7] H.K. Park, M.Y. Ha, H.S. Yoon, Y.G. Park, Changmin Son, A numerical study on natural convection in an inclined square enclosure with a circular cylinder, Int. J. Heat Mass Transfer 66 (2013) 295–314.

[19] H.S. Yoon, Y.G. Park, J.H. Jung, Natural convection in a square enclosure with differentially heated two horizontal cylinders, Numerical Heat Transfer 65 (2014) 302–326. [20] Y.G. Park, M.Y. Ha, C. Choi, J. Park, Natural convection in a square enclosure with two inner circular cylinders positioned at different vertical locations, Int. J. Heat Mass Transfer 77 (2014) 501–518. [21] Y.G. Park, H.S. Yoon, M.Y. Ha, Natural convection in square enclosure with hot and cold cylinders at different vertical locations, Int. J. Heat Mass Transfer 55 (2012) 7911–7925. [23] E.M. Sparrow, J.E. Niethammer, Effect of vertical separation distance and cylinder-to-cylinder temperature imbalance on natural convection for a pair of horizontal cylinders, J. Heat Mass Transfer 103 (1981) 638–644. [29] J. Kim, P. Moin, Application of a fractional step method to incompressible Navier-Stokes equations, J. Comp. Phys. 59 (1985) 308–323. [30] Y. Zang, R.L. Street, J.R. Koseff, A non-staggered grid, fractional step method for time-dependent incompressible Navier-Stokes equations in curvilinear coordinates, J. Comp. Phys. 114 (1994) 18–33.