Effect of spacing on laminar natural convection flow and heat transfer from two spheres in vertical arrangement

International Journal of Heat and Mass Transfer 134 (2019) 852–865

Contents lists available at ScienceDirect

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

Effect of spacing on laminar natural convection flow and heat transfer from two spheres in vertical arrangement Jian Zhang a, Qi Zhen b, Jie Liu a,⇑, Wen-Qiang Lu a a b

School of Engineering Science, University of Chinese Academy of Sciences, Beijing 100049, China College of Energy and Traffic Engineering, Inner Mongolia Agricultural University, Hohhot 010018, China

a r t i c l e

i n f o

Article history: Received 25 October 2018 Received in revised form 17 December 2018 Accepted 16 January 2019

Keywords: Numerical investigation Natural convection Two vertically arranged spheres Sphere spacing Heat transfer correlations

a b s t r a c t Laminar natural convection heat transfer from two vertically arranged spheres has been investigated numerically in air (Pr ¼ 0:72) for large Rayleigh number (Ra ¼ 105 ). The effect of sphere spacing (S) on fluid flow and heat transfer is comprehensively analyzed within the ratio of spacing to radius (S=R) ranging from 4 to 48, i.e., 2 6 S=D 6 24. Extensive results are obtained for complete descriptions on natural convection of two-sphere system, such as the streamline patterns and isotherm contours, dimensionless velocities and temperature profiles, drag coefficient, local and average Nusselt numbers. Comparisons with the isothermal bisphere have also been reported with the previous studies. Due to the interaction of preheating effect and flow mixing, the lower sphere has a great influence on the fluid flow and heat transfer around the upper sphere with the change of spacing between the two spheres. The variation of dimensionless temperature with the dimensionless distance between the two spheres has been also concluded to explain the impact of preheating effect. Finally, correlations have been developed for accurate prediction of natural convection heat transfer for each single sphere and the whole two-sphere system, which can be applied to relevant engineering calculations, and provide values for academic research. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Natural convection is a flow motion not driven by external forces (except for the gravity) but caused by the inhomogeneity of the temperature or concentration fields within the fluid. This inhomogeneity leads to the appearance of density gradient, and the resulting buoyancy becomes the driving force of fluid motion [1]. Heat transfer process through natural convection is prevalent in industry because of its inherent advantages, such as economy, adequate safety and low noise, etc. As one of the most fundamental structures, the spherical geometry often occurs in all kinds of industrial or engineering applications. The phenomenon of natural convection heat transfer from spherical body is always encountered in the following fields of heat transfer [2–5], such as spherical electric lights, fuel droplets, fixed beds, fluidized beds, chemical reactors, pebble bed nuclear reactor, foodstuffs treatment and so on. Undoubtedly, the heat transfer of hot spheres (or particles) to fluids has become a critical research, which is of great research value. Many scholars have made detailed studies on the natural convection heat transfer of a single sphere in the previous studies, ⇑ Corresponding author. E-mail address: [email protected] (J. Liu). https://doi.org/10.1016/j.ijheatmasstransfer.2019.01.065 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

from the aspects of theory, experiment and simulation [6–12]. Some of them have also given empirical expressions of natural convection to facilitate the prediction of heat transfer. In practice, however, heat transfer from a single sphere is not the most common situation. The heat transfer phenomenon of double spheres or multi-spheres is more universal and representative, because it considers the interaction of adjacent structures in a system. Unfortunately, there are very few studies focusing on this kind of problem due to its complexity. On the other hand, for spheres with different spacing, the characteristics of natural convection heat transfer and fluid flow are also seldom studied. In early research, Chamberlain et al. [13] conducted experimental and theoretical study on the natural convection heat transfer from an isothermal bisphere (two touched spheres) in air. Based on the R-H method, they presented a correlation with an average deviation of 3% for predicting the heat transfer rate from a bisphere for 10 6 Ra 6 107 . According to the linear superposition of the theoretical diffusive limit (Ra ! 0) and thin boundary layer asymptomatic limit (104 < Ra < 108 ), natural convection from complex isothermal three-dimensional bodies (including the bisphere) has been theoretically analyzed by Yovanovich [14]. On the basis of one empirical parameter (characteristic body length), a semiempirical correlation with the accuracy equal to the complex

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865

853

Nomenclature CD C DF C DP D D1 g Gr h k L Nu n p Pr R Ra s S T Ts Tm DT u; v

total drag coefficient (dimensionless) friction drag coefficient (dimensionless) pressure drag coefficient (dimensionless) diameter of the sphere (m) width of the computational domain (m) gravitational acceleration (m=s2 ) Grashof number (dimensionless) convective heat transfer coefficient (W=m2  K) thermal conductivity (W=m  K) radial distance from the surface of the sphere (m) average Nusselt number (dimensionless) power-law index (dimensionless) gauge pressure (Pa) Prandtl number (dimensionless) radius of the sphere (m) Rayleigh number (dimensionless) surface area of the sphere (m2 ) spacing between the center of the two spheres (m) temperature of the fluid (K) surface temperature of the sphere (K) film temperature of the fluid (K) temperature difference, (¼ T s  T 0 , K) x and y components of the velocity (m=s)

equations of Raithby-Hollands [15] and Hassani-Hollands [16] is also proposed for the rough prediction of total Nusselt number. Then they experimentally investigated the natural convection heat transfer from three-dimensional bodies of arbitrary shape (including the bisphere) in air (Pr  0:71) for Rayleigh number ranging from 10 to 108 [17]. By simplifying the Raithby-Hollands method, they developed a general expression for predicting the heat transfer rate with an average rms difference of 3.5%. In recent years, study on the natural convection from two sideby-side or tandem spheres has begun to emerge. Leweke et al. [18] studied the interactions of two spheres placed side by side in a uniform stream at low Reynolds numbers by experiments and direct numerical simulation (DNS). Flow characteristics of two side-byside spheres immersed in a uniform flow at Re ¼ 5  103 have been studied experimentally within the spacing range of 1:0 6 S=D 6 2:50 [19]. For two tandem spheres as an array immersed in power-law fluids, Mukherjee et al. [20] explored the effect of the power-law index (0:4 6 n 6 1:6) and inter-particle distance (2 6 S=D 6 8) on local velocity and temperature fields. By employing the three-dimensional immersed boundary method (IBM), Musong et al. [21] investigated numerically the thermal interactions between spherical particles in a viscous Newtonian fluid, natural convection heat transfer from bispheres in tandem and in line has been considered in air for small Grashof number (Gr ¼ 100). The heat transfer correlations obtained from the relevant studies are listed in Table 1. Obviously, it has become more and more important to study the interactions of natural convection heat transfer and flow structure from two-sphere system, especially a wide range of inter-sphere spacing and large Grashof or Rayleigh number (Ra ¼ Gr  Pr). However, the previous studies have not paid much attention on this issue. As a result, this work has investigated thoroughly the natural convection heat transfer and fluid flow from two vertically arranged spheres in air. The effect of inter-sphere spacing and large Rayleigh number on the streamline and isotherm contours, velocity and temperature fields, drag coefficient, as well as Nusselt numbers have been studied for 2 6 S=D 6 24 and Ra ¼ 105 .

U; V Uc Vr Vh x; y X; Y

x and y components of the velocity (dimensionless) reference velocity (m=s) radial velocity (m=s) tangential velocity (m=s) Cartesian coordinates (m) Cartesian coordinates (dimensionless)

Greek symbols a thermal diffusive coefficient (m2 =s) b coefficient of thermal expansion (K1 ) h circumferential angle (degree) # kinematic viscosity (m2 =s) q density of the fluid (kg=m3 ) l fluid viscosity (Pa  s) d grid size on the surface of the sphere (m) Superscript and subscripts  dimensionless quantities 0 ambient condition h local value L lower sphere U upper sphere W the whole collection

The obtained results have been analyzed in detail, which are in good agreement. Meanwhile, steady-state heat transfer correlations for predicting respective heat transfer of single sphere and overall heat transfer of two-sphere system are presented by data fitting. 2. Problem formulation In this work, two vertically arranged spheres of diameter D maintained a uniform surface temperature of T s , are immersed in an extensive and quiescent air of T 0 . Obviously, the fluid near the surface of the sphere is heated because of T s > T 0 , hence the density difference caused by the temperature gradient will induce a steady upward buoyant flow within the region, accompanying by corresponding convective heat transfer. As shown in Fig. 1, the considered structure is placed in a channel, the distance between the center of two vertically arranged spheres is S, and the origin of the Cartesian coordinate system is set at the midpoint of the centerline between the two spheres. The computational model based on Bejan et al. [22] is employed in this calculation, which has been successfully applied in two attached cylinders [23,24] and hemisphere [25,26]. For Ra ¼ 105 , it is perfectly reasonable to consider the flow as laminar, steady-state and axisymmetric [1]. Therefore, the traditional three-dimensional model of the sphere can be simplified to a two-dimensional situation. Meanwhile, the thermophysical properties are assumed to be temperature independent as the temperature difference is small (DT < 5 K). As for the density item in y-direction momentum equation, it is computed by the film temperature, T m ¼ ðT s þ T 0 Þ=2. Also, the Boussinesq approximation, q0  q ¼ qbðT  T 0 Þ, is used to evaluate the variation of density with temperature, where b is the thermal expansion coefficient acquired by b ¼ 1=qð@ q=@T Þp . For small temperature difference case, b can be simply written as b ¼ 1=T m [1]. Without considering the viscous dissipation and radiation (due to the low velocity and small temperature difference in natural convection), the governing equations of the present flow are simplified as follows:

854

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865

Table 1 Summary of natural convection heat-transfer correlations for two-sphere system. Authors

Correlations

Chamberlain et al. [13]

Nu ¼



Restrictions  1 4:8

1:38 þ 0:378Ra4

 1=4:8 1 4:8 þ 0:104Ra3 

Bisphere, air, 0 < Ra < 107 , experimental

Yovanovich [14]

Nu ¼ 3:470 þ 0:510Ra1=4

Bisphere, air, 0 < Ra < 108 , experimental

Raithby and Hollands [15]

NuT ¼ 0:745C l Ra1=4

Bisphere, air, 1 < Ra < 107 , experimental

Cl ¼ 

0:671 9

1þð0:492=PrÞ16

4=9

  1:02 1=1:02 Nul ¼ 1:391:02 þ NuT

Hassani and Hollands [17]

Nut ¼ 0:104Ra1=3 h i1=15 15 Nu ¼ Nu15 l þ Nut  1=1:07  1 1:07 þ 3:4751:07 Nu ¼ 0:114Ra4

Bisphere, air, 10 < Ra < 108 , experimental

Schouveiler et al. [18] Pinar et al. [19]

None None

Two side-by-side spheres, water, Re ¼ 300, 1 6 S=D 6 3, numerical

Mukherjee et al. [20]

Nu ¼ 1:725X0:682 ð1 þ ð1=ðS=DÞÞÞ0:778 (X: user-defined parameter) None

Two side-by-side spheres, power-law fluids, 102 6 Ra 6 105 , 2 6 S=D 6 8, numerical

Musong et al. [21]

Two side-by-side spheres, water, Re ¼ 5  103 , 1 6 S=D 6 2:5, experimental

Two side-by-side spheres, air, Gr ¼ 100, 1 < S=D < 10, numerical

Fig. 1. Schematic representation of numerical model and computational domain.

Continuity equation:

@u @ v þ ¼0 @x @x

ð1Þ

Momentum equations:

x-direction :

u

@u @u 1 þv ¼  rpx þ #r2 u @x @x q

ð2Þ

y-direction :

u

@v @v 1 þv ¼  rpy þ #r2 v þ gbðT  T 0 Þ @x @x q

ð3Þ

Energy equation:

u

@T @T þv ¼ ar2 T @x @x

ð4Þ

where r2 ¼ @ 2 =@x2 þ @ 2 =@y2 , and p is the gauge pressure, u, v are the velocities in x and y directions, #, a are the kinematic viscosity and thermal diffusivity, respectively. In addition, the definitions of all variables appearing in Eqs. (1)– (4) are non-dimensionalized by the following formulas:

X¼

x y T  T0 u v  pR ; Y ¼ ; T ¼ ;U ¼ ;V ¼ ;p ¼ R R Uc Ts  T0 Uc lU c

ð5Þ

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where U c is the reference velocity defined by U c ¼ gbDTR, l is the fluid viscosity. Some dimensionless numbers are involved in studying the problem of natural convection heat transfer of the two spheres, which are used to analyze the flow and heat transfer. For the

855

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865 Table 2 Total number of grids varying with different spacing ratios (S=D). S=D

2

4

8

16

24

Total number of grids

512,100

700,906

1,108,662

1,904,340

1,925,746

convenience of understanding, their definitions are introduced as follows: Rayleigh number (Ra):

Ra ¼ Gr  Pr ¼

gbDTR3 #a

ð6Þ

Prandtl number (Pr):

Pr ¼

#

a

¼

Cp l k

ð7Þ

where C p is the heat capacity at constant pressure, k is the thermal conductivity. Local Nusselt number (Nuh ):

Nuh ¼



hh R @T ¼ k @ns surface

ð8Þ

where ns is the unit normal vector perpendicular to the surface of the sphere. Average Nusselt number (Nu):

Nu ¼

Z  hR 1 ¼ Nuh ds k s s

ð9Þ

Drag coefficient (C D ): It consists of two parts, namely the friction drag coefficient (C DF ) and the pressure drag coefficient (C DP ), which are caused by the shear force and normal force, respectively. Their respective definitions are as follows:

C D ¼ C DF þ C DP C DF

ð10Þ

Z  1 ¼ F D;f qU 2c pR2 ¼ C f  ns ds 2 s

ð11Þ

where F D;f is the frictional drag, C f is the local friction drag coefficient defined by:

,



1 2 @U h qU 2c ¼ pffiffiffiffiffiffi 2 Gr @ns surface

ð12Þ

Z  1 ¼ F D;p qU 2c pR2 ¼ C p  ny ds 2 s

ð13Þ

C f ¼ F D;f

C DP



1 qU 2c 2

have been carried out with half of the computational domain, namely, the right half region of X P 0. As shown in Fig. 1, the corresponding boundary conditions are given out, and their brief definitions are listed as follows: On the sphere surface: no-slip and isothermal boundary conditions are prescribed, i.e.,

u ¼ v ¼ 0; T ¼ T s :

where F D;p is the pressure drag given by F D;p ¼ Ps  P 0 , ny is the unit vector perpendicular to the surface in the y direction, C p is the local pressure drag coefficient defined by:

C p ¼ F D;p

Fig. 2. Diagrammatic sketch of computational grids around the two-sphere system.

ð14Þ

3. Numerical methodology According to the above descriptions, numerical investigations on steady natural convection from two vertically arranged spheres

ð15Þ

The first inlet: the fluid flows into the computational domain from the bottom, which is set as the ‘‘inlet”. The normal stress, vertical velocity gradient and horizontal velocity on this surface are all zero (to keep the flow in a perfectly pure natural convection, without interference from external forces except gravity), and the fluid temperature is set to the ambient temperature,

p ¼ 0; @ v =@y ¼ 0; u ¼ 0; T ¼ T 0 :

ð16Þ

The second inlet: to avoid the chimney effect (cause artificial fluid acceleration), both sides of the outlet section (region C) are enforced as the ‘‘inlet”. The difference is that the horizontal velocity gradient and vertical velocity are zero,

p ¼ 0; @u=@x ¼ 0; v ¼ 0; T ¼ T 0 :

ð17Þ

Table 3 Domain independence tests at Gr ¼ 107 . Domain length Inlet section

Outlet section

D1

C D;p

CD

Nu

3D 6D

3D 6D

3D 6D

0.0145 0.0148

0.0292 0.0293

24.904 24.868

856

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865

Table 4 Grid independence tests at Gr ¼ 107 . Grid

Total number of cells

Smallest cell size d=D

C D;p

CD

Nu

Grid 1 Grid 2

980,202 2,791,076

0.001 0.0005

0.0145 0.0145

0.0292 0.0293

24.904 24.878

Fig. 3. Comparisons of average Nusselt number for biosphere with the previous studies.

Outlet: the fully flowing fluid comes out from the top of the computational domain, which is set as the ‘‘outlet”, and the backflow temperature is also set to the ambient temperature,

p ¼ 0; @ v =@y ¼ 0; u ¼ 0; T ¼ T 0 :

ð18Þ

At the centerline: due to the axisymmetric flow for the considered Rayleigh number, the centerline is set as the ‘‘axisymmetry”,

@ v =@x ¼ 0; u ¼ 0; @T=@x ¼ 0:

ð19Þ

At the sides: to reduce the size of computational domain, the two sides of the inlet section (region A) plus the sphere section (region B) are defined as the ‘‘symmetry”,

@T=@x ¼ @ v =@x ¼ 0; u ¼ 0:

ð20Þ

As for the steady-state numerical solution, the governing equations mentioned in the previous section combined with these boundary conditions have been solved by using the ANSYS Fluent (version 15.0), which are iteratively solved by the control volume formulation and discretized by the finite difference method (FDM), respectively. The SIMPLE scheme is applied to dealing with the pressure-velocity coupling. For convective terms in momentum and energy equations, the QUICK scheme has been chosen. Under the consistent accuracy, the Least Squares Cell Based method is employed in gradient discretization to save the computational efforts (Due to the large amount of meshes, especially when the spacing is large). Moreover, the convergence criterions for continuity, momentum and energy equations are of the order of 106 in order to obtain precise results. To learn more about the numerical simulation procedure, it can be found in our latest studies [25,26]. In general, numerical simulation on natural convection needs a huge space to avoid the effect caused by the boundary. As noted in the study of Prhashanna and Chhabra [10], the ratio of D1 =D is at least 80, the boundary has no influence on the fluid flow and heat transfer. As a result, in the previous studies, the computational

width is always more than 80 times of the diameter for spheres [9,10,27], 400 times for a cone [28] and even 600 times for cylinders [29], which leads to enormous computational efforts. However, an innovative model has been proposed by Bejan et al. [22], which is used to conduct simulation on the natural convection from horizontal cylinders in a fixed volume. With the help of this model, the modified one is applied in the current calculation of natural convection from two spheres, as shown in Fig. 1. The advantage of this model is that it can effectively reduce the computational domain (at least 10 times) and save the numerical costs while ensuring the accuracy of the calculation results, thus improving the computational efficiency. Table 2 shows the total number of grids varying with different spacing ratio (S=D) in the present study. It can be seen that the total number of grids for S=D ¼ 24 is almost four times than that of S=D ¼ 2, although the present model can greatly reduce the size of computational domain. In our previous works [25,26], the accuracy and reliability of this model have been fully verified, except for some different boundary conditions. To find appropriate numerical parameters, the domain dependence tests for single sphere have been done for high Grashof number of 107 (i.e., Ra ¼ 7:2  106 ). According to our previous works, the computational domain is reasonable enough to choose 3D  3D  3D and 6D  6D  6D. Table 3 shows that the results of C DP , C D and Nu for two situations have hardly changed, although the computational regions for the inlet section, outlet section and D1 have doubled. As a result, with this high conformance, selecting the numerical parameters of 3D  3D  3D is reliable and accurate, which is much smaller than the traditional model. On the other hand, the computational grids adopted in this work consists of two part, namely, the highly fine mesh near the sphere surface (to capture the sharp variations of velocity and temperature) and the relatively coarse mesh far away from the sphere (to economize the computational costs). As shown in Fig. 2, the mesh of part I is finely meshed by the non-uniform quadrilateral grids with a small number of triangular grids, while the mesh of part II is generated by the structural quadrilateral grids with an expansion ratio of 1.02 in the x and y directions. It is known that the preferable mesh is required to deal with the dramatic change in velocity and temperature gradients at large Grashof or Rayleigh number (Ra ¼ Pr  Gr). Thus, the grid independence tests for single sphere have also been carried out for high Grashof number of 107 . As demonstrated in Table 4, the G1 has a total of 980,202 grids with the smallest grid size d=D ¼ 0:001, while the G2 has 2,791,076 grids with d=D ¼ 0:0005. Although the total grids number of G2 is about three times that of G1, the values of C DP and C D remain almost unchanged, the value change of Nu is also just about 0.34%. Consequently, the choice of G1 with d=D ¼ 0:001 is considered to be fine enough in this study. In a word, it is perfectly reasonable to choose the computational domain width of 3D and the mesh with grid density of 0.001 as the numerical parameters. What’s more, the numerical solution applied in this work should be further validated to ensure its reliability and accuracy. Comparisons with the previous works have been made in the next section.

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865

857

Fig. 4. Isotherm contours (left half) and streamline (right half) around the two-sphere system for S=D varied from 2 to 24 at Ra ¼ 105 .

4. Results and discussions

spacing between the two separate spheres on the fluid flow and heat transfer has been investigated for 2 6 S=D 6 24 and Rayleigh

Numerical simulation on the natural convection from two vertically arranged spheres has been conducted in air. The effect of

number of 105 . However, before showing the obtained results, the accuracy and reliability of the new results should be further veri-

858

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865

4.1. Validation of results So far, natural convection heat transfer of two vertically arranged spheres in close contact has been studied by Chamberlain et al. [13], Hassani and Hollands [17]. In this validation, the average Nusselt number of an isothermal bisphere has been obtained for 10 6 Ra 6 107 in air by using the new numerical methodology. The errors of average Nusselt number are compared to verify whether it is suitable for studying the natural convection heat transfer of bisphere. As shown in Fig. 3, the error of numerical results within 102 to 106 is in good agreement with the results of previous studies, and the maximum error is no more than 2%. At Ra ¼ 10, the numerical results are closer to the theoretical results of Hassani and Hollands [17] than those of Chamberlain et al. [13]. At Ra ¼ 107 , the results are between the results of above studies, and the maximum error is still less than 3.5%. Accordingly, the results are completely acceptable and prove the applicability of our calculation model in a wide range of Rayleigh numbers for the natural convection heat transfer of two-sphere system. 4.2. Isotherm contours and streamline patterns

Fig. 5. Dimensionless velocity curves of x and y directions in the range of 105—172:5 .

To demonstrate qualitatively the effect of different spacing on the temperature and velocity distributions of two side-by-side spheres, the isotherm contours (left half) and streamline patterns (right half) for 2 6 S=D 6 24 have been shown in Fig. 4(a)–(g). It can be seen that even though the ratio of distance to diameter is 24 times the sphere diameter, the lower sphere still has a slight effect on the fluid flow and heat transfer of the upper sphere. The convective flow that rises from the lower sphere is mixing with the upper one, and forming into the boundary layer. As heated by the lower sphere, the fluid gets a certain velocity when passing the upper sphere, which makes the fluid flow and heat transfer around the upper sphere enhanced by the mixing effect of convection. It can be seen from Fig. 4(a)–(g) that when the distance between the two spheres is small, the streamlines near the sphere surface are deflected to a large extent due to the effect of entrainment between the two spheres. However, with the increase of distance, the deflection of the streamlines becomes less obvious, and the velocity boundary layer near the sphere surface is gradually thickened. It shows that the velocity change of the plume increases the strength of the fluid flow when the distance between the two spheres is small, which leads to the enhancement of heat transfer of the two spheres. Further observation shows that the movement of the fluid along the lower sphere is very similar to that of a single sphere. However, a recirculation zone is gradually formed at the top of the upper sphere when the ratio of S=D P 20. This is due to the fact that the increase of distance produces the characteristics of mixed flow when the fluid flows through the upper sphere. Fig. 5 shows the dimensionless velocity curves of x and y directions at S=D ¼ 24, which can be seen that the velocities (Fig. 5(a) and (b)) 

Fig. 6. Setting and schematic of related parameters in data processing (angle, radial and tangential velocities).

fied. Therefore, comparisons of theoretical and experimental results between the present results and the available studies in the literature have been established.



have reversed in the range of 142:5  180 when the dimensionless radial distance (¼ L=R) is less than 1.41 (R 6 1:41). Due to the heating by the sphere, the viscosity of fluid on the sphere surface is largest, and rapidly decreases along the radial direction of the sphere. This is equivalent to the formation of a stable layer of viscous fluid enclosing the sphere at a very close distance. For two vertically arranged spheres, the viscous thin layer formed by the fluid consists of the adjacent region of the spheres and the connection region between the two spheres. On the other hand, since the flow direction of the fluid is vertical, which is equivalent to enlarging the effective domain of the upper sphere immersed in the fluid, thus increasing the local Nusselt number at the top of the sphere, and leading to the formation of recircula-

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865

859

Fig. 7. Radial velocity distribution along various directions for S=D ¼ 4 and 20.

tion zone over the top (Fig. 4(g)). As seen from Fig. 4, with the increase of the distance between the two spheres, the increasing volume of the viscous thin layer makes the effective size of the upper sphere bigger (which can be seen in Fig. 9), and the point of flow separation around the upper sphere will move upstream, which brings about the larger recirculation zone. 4.3. Velocity and temperature distributions As shown in Figs. 7–9, the dimensionless radial velocity (¼ V r =U c ) and tangential velocity (¼ V h =U c ), temperature (¼ ðT  T 0 Þ=ðT s  T 0 Þ) of each sphere have been investigated quantitatively along various angles for small ratio of S=D ¼ 4 and large ratio of S=D ¼ 20. The influence of sphere spacing on the distributions of above parameters are analyzed in detail, and the setting of related parameters in data processing are sketched in Fig. 6. Fig. 7 shows the variation of radial velocity along the dimensionless radial distance (R ¼ L=R) at various angles. In general, it can be divided into two parts by h ¼ 90 . After R > 2, the values of h ¼ 90 are almost zero (a little bit of shock is caused by the interactions), however, the values of h ¼ 0—90 are negative while the values of h ¼ 90—180 are positive, which can be attributed to the curvature effect. For h < 90 , the negative value means that the

fluid flows toward the sphere and the domain can be regarded as the inlet (this is consistent with the setting of boundary conditions in computational domain). On the contrary, for h > 90 (especially at large angles), the flow direction has changed and the fluid flows away from the sphere (so that the velocity is positive), and these regions can be regarded as the outlet. Generally, the value of small angle first decreases and then stabilizes at a level, except for h ¼ 0 at S=D ¼ 4 due to the influence of the lower sphere (Fig. 7(a)). For 90 < h < 135 , it increases first and then decreases. When S=D ¼ 4, the values of the lower sphere have a rapid drop at h ¼ 135 and 180 for R P 6 (Fig. 7(b)) because of the flow interactions within the two spheres. For h ¼ 180 , the fluid continues to flow away from the sphere, so its value rises slowly and eventually reach to a level. As for the variation of tangential velocity seen in Fig. 8, the values of h ¼ 0 and 180 are basically keeping at the same level for both spheres. The general trend of the other curves is to increase rapidly and then to decline slowly, except for h ¼ 135 of the lower sphere at small spacing S=D ¼ 4 (Fig. 8(b)), resulting in the influence of the upper fluid. Meanwhile, it takes longer for the upper sphere (R ¼ 5) to reach stability than for the lower one (at R ¼ 4) under the interactions of upcoming flow. The maximum ratio of the two spheres (V h;U =V h;L ) is about 1.13 for S=D ¼ 4, and

860

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865

Fig. 8. Tangential velocity distribution along various directions for S=D ¼ 4 and 20.

3.67 for S=D ¼ 20, respectively, which shows the effect of sphere spacing on its distribution. It can be found that with the increase of spacing, the velocity gradient of each sphere will increase relatively because of the thinning boundary layer. In Fig. 9, all the variations of temperature decrease rapidly with the increasing radial distance, which are close to zero at R ¼ 2, except for h ¼ 180 for both spheres and h ¼ 0 for the upper sphere. For the upper sphere, the value of h ¼ 0 keeps going up (Fig. 9(a)) due to the heating effect by the lower one at small spacing (S=D ¼ 4), while the heating effect of the lower sphere is stable at large spacing (S=D ¼ 20), and its value almost remains at a horizontal level (Fig. 9(a)). The value of h ¼ 180 is always larger because the heat is constantly flowing upward during natural convection. On the other hand, the descending rate of the upper sphere is smaller than that of the lower one, which causes a longer time to reach zero in the whole process (about R ¼ 2 for the lower sphere and R ¼ 3 for the upper sphere). 4.4. Analysis of local Nusselt number Based on the analysis of local Nusselt number, a deeper explanation of the spacing effect between the two spheres on the fluid flow and heat transfer is obtained for 2 6 S=D 6 24. Fig. 10(a)

shows the distribution of local Nusselt number at different positions on the surface of the lower sphere. On the whole, the value of the lower sphere reaches the maximum at the front stagnation point. As the fluid moves along the sphere surface, the curve decreases monotonously, and attains the minimum value at rear stagnation point. The curve trend is similar to that of the isothermal single sphere, which is caused by the gradual thickening of the temperature layer along the sphere surface and the constant decay of the temperature gradient. However, the maximum local Nusselt number of the lower sphere is positively dependent on the distance between the two spheres, which is usually determined by the temperature gradient on the sphere surface (usually affected by the thickness of the temperature field and the heat transfer rate of the fluid). From the variation of temperature field in the previous section, it can be seen that the temperature field of the lower sphere bottom does not become thinner with the increase of distance between the two spheres. Nevertheless, the larger the distance is, the more fluid can be touched with the sphere for heat transfer, which makes the heat transfer capacity of the lower sphere stronger. For h ¼ 135 , the fluid begins to separate around the lower sphere. After this point, the temperature gradient is more and more affected by the thickness of the temperature boundary layer, and the influence of the fluid heat transfer

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865

861

Fig. 9. Temperature distribution along various directions for S=D ¼ 4 and 20.

rate is minimal. When the angle is greater than 135 , the curves with different spacing eventually fall into one curve because the boundary layer becomes thicker as the angle increases. The distribution of local Nusselt number on the surface of the upper sphere is shown in Fig. 10(b). As mentioned above, the heat transfer rate of the upper sphere is affected by the lower one, and its heat transfer capacity is greatly affected by different distances between the two spheres. For S=D ¼ 2, the curve trend is different from the others. The curve keeps rising from h ¼ 0 to h ¼ 20 , and reaches a maximum in the vicinity of h ¼ 40 , then slowly descends to the first confluence point at h  120 . The curve drops rapidly for h ¼ 120—180 , and ends at h ¼ 180 with getting a minimum. In the process, there is the second confluence point at h  160 . As we can see, when the distance between the two spheres is small, the fluid contacting with the lower part of the upper sphere is preheated by the lower one, which makes the position of the forward stagnation point move to h ¼ 40 . On the other hand, as the spacing increases, more cold fluid is carried to the vicinity of the sphere from the outside, accompanying with heat exchange between cold and hot fluids. This makes the preheating of the lower sphere count for little, and the heat transfer enhancement of the upper sphere is due to the contribution of the resulting mixed flow. Therefore, the front stagnation point is also the position of the largest local Nusselt number when the value of S=D is

large. With the increase of recirculation zone, the flow separation point of the upper sphere shifts to the upstream direction, which makes the point of the minimum value of local Nusselt number move forward gradually, and the increasing mixed motion also increases the minimum value of local Nusselt number. 4.5. Drag coefficient Due to the normal force and tangential force applied to the fluid (resulting from the temperature difference between the two spheres), the fluid produces a velocity gradient, and the resultant force is usually expressed by the drag coefficient. In the previous studies, some of them have studied the heat transfer coefficient of various shapes. However, there are few literatures to analyze the influence of the fluid resistance for the system of two spheres. Fig. 11 shows the dependence of average drag coefficient on S=D for the two spheres with various spacing. The average drag coefficient of the upper sphere increases with the increasing distance between the two spheres, but the curve of the lower sphere remains at the same level. It proves that the increasing mixing motion near the upper sphere causes the fluid to produce a sharp increase in velocity gradient, resulting in greater resistance than the lower sphere. Fig. 12 shows the ratio of the pressure drag coefficient to the total drag coefficient of the two spheres, which is to determine

862

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865

Fig. 12. Variation of the ratio C DP =C D versus S=D for each sphere in the system.

Fig. 10. Variation of local Nusselt number along the surface of each sphere for S=D varied from 2 to 24.

Fig. 13. Variation of dimensionless temperature with the distance at the centerline.

Fig. 11. Variation of average drag coefficient versus S=D for each sphere in the system.

Fig. 14. Correlations of average Nusselt number as a function of S=D for each sphere in the system.

863

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865 Table 5 Comparisons between numerical results and formula fitting results. S=D

Numerical results

Formula fitting results

Error

Lower sphere 2 9.816 4 9.854 16 10.089 24 10.246

9.815 9.858 10.102 10.250

0.02% 0.04% 0.13% 0.04%

Upper sphere 2 8.305 4 10.158 16 14.066 24 14.961

8.432 10.017 14.139 14.818

1.51% 1.41% 0.52% 0.96%

the centerline. It can be concluded that the influence of the lower sphere on the flow and heat transfer of the upper sphere becomes smaller with the increasing spacing. When S=D > 12, the preheating effect has been softened. And there is just a weak impact while S=D ¼ 24. Fig. 14 shows the variation of average Nusselt number for the two spheres with different spacing. The exact correlation of average Nusselt number as a function of S=D is obtained for the large Rayleigh number (Ra ¼ 105 ). It can be seen from the diagram that when 2 6 S=D < 4, the average Nusselt number of the lower sphere is higher than that of the upper one. This is because the preheating effect under small gap reduces the heat transfer rate of the upper sphere (this is equivalent to the reduction of effective temperature difference). However, as the distance between the two spheres increases up to 4D, the average Nusselt number of the upper sphere turns to reverse, and constantly widen the gap. Due to the increasing distance between the two spheres and the virtual mixed flow around the sphere, these factors balance the adverse effect of the preheating by the lower sphere, and thus improve the heat transfer capacity of the upper sphere. According to the numerical results, the respective correlations of average Nusselt number are obtained as follows: For the lower sphere:

NuL ¼ 12:1417  2:3712  0:9906ðS=DÞ

ð21Þ

For the upper sphere:

NuU ¼ 15:1716  8:8108  0:8746ðS=DÞ

Fig. 15. Correlation of average Nusselt number versus S=D for the whole collection.

the role of friction and pressure resistance in the process of increasing spacing. The ratio of the lower sphere keeps a slow decline for the whole range of S=D. Meanwhile, as the distance between the two spheres increases, more cold fluid swarms into the gap region from both sides, which slightly improves the heat transfer ability of the fluid to the lower sphere, increases the velocity gradient of the fluid and slowly enlarges the friction resistance of the fluid. For 2 6 S=D 6 4, the ratio of the upper sphere decreases gradually, which indicates that the preheating effect of the lower sphere reduces the velocity gradient of the fluid near the sphere, and increases the corresponding friction resistance. This is also consistent with the fact that the average Nusselt number of the upper sphere is lower than that of the lower one in the range of 2 6 S=D 6 4. When the spacing between the two spheres exceeds 4D (S=D > 4), the preheating effect is weakened, and the mixing flow enhances the pressure gradient of the fluid near the upper sphere, which accelerates the velocity of the fluid, so the ratio of C DP =C D increases with the increasing spacing. As a result, the difference in C DP =C D of the two spheres is gradually widening.

ð22Þ

The determination coefficients of Eqs. (21) and (22) are 0.9362 and 0.99786, respectively. Table 5 shows the numerical results and the calculated results by the fitting formula, respectively, and the relative errors between them. As we can see, the average error of the lower sphere is about 0.088%, and the maximum error is less than 0.173%. For the upper sphere, the average error of is about 0.785%, and the maximum error is 1.51%. This shows that the formula is in good agreement with the numerical results and has a high accuracy. On the other hand, as shown in Fig. 15, we add up the average Nusselt number of the upper and lower spheres over the range of 2 6 S=D 6 24, and obtain the correlation of the whole two-sphere system by fitting the data:

NuWhole ¼ 25:5812  9:40  0:8839ðS=DÞ The

determination

coefficient

ð23Þ of

the

above

equation,

R2 ¼ 0:99805, which illustrates that the existing results have enough precision. Moreover, all the formulas obtained above can offer accurate prediction of natural convection heat transfer from a single sphere or the whole in two-sphere system under the same or similar conditions, which has certain industrial value and academic significance.

5. Conclusions

4.6. Correlations of average Nusselt number

In this paper, we investigated numerically the natural convection heat transfer of the two spheres in vertical arrangement for 2 6 S=D 6 24 (i.e., 4 6 S=R 6 48) at large Rayleigh number

In general, the effect of spacing on the heat transfer of the two spheres can be analyzed in detail by the local Nusselt number of each sphere. Meanwhile, the average Nusselt number of a single sphere or the whole is often encountered in the design and calculation of particle heat exchanger and so on. As shown in Fig. 13, The preheating effect of the lower sphere on the upper one has been revealed through the variation of dimensionless temperature with the distance between the surface of the two spheres (D 2 ½0; 1 ) at

(Ra ¼ 105 ). The distributions of temperature and velocity are analyzed in detail. Meanwhile, the local Nusselt number distribution on the surface of each sphere is given out, and the variation of their respective average Nusselt numbers with different spacing is summarized. The preheating effect of the lower sphere on the upper one has been explored through the normalized temperature variation curve at the centerline. By analyzing the average drag coefficient of each sphere, the characteristics of fluid flow near the

864

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865

two spheres are also presented. The main conclusions are summarized as follows: (1) For two spheres with spacing, the streamline around the surface has a large degree of deflection due to the entrainment of the two spheres at small S=D. For the lower sphere, its flow characteristic near the surface is very similar to that of a single isothermal sphere. As the spacing increases to 24 times the diameter, the recirculation zone will occur over the upper sphere for 142:5 6 h 6 180 when R 6 1:41 because of the mixing effect. And the flow separation point around the upper sphere moves upstream with the increasing spacing between the two spheres. For S=D ¼ 4, the radial velocity of the lower sphere has a rapid drop at h ¼ 135 and 180 for R ¼ 6 due to the flow interactions. Besides, for the upper sphere, it takes longer to reach a stable level in the distributions of radial and tangential velocities, and temperature. (2) The convective heat transfer coefficient of the upper sphere depends on the two competing factors: the preheating effect and the flow mixing from the lower one, which are seriously affected by the value of S=D. When S=D < 4, the effect of preheating on the heat transfer of the upper sphere is greater than that of the mixed flow, and the average Nusselt number of the upper sphere is lower than that of the lower one. When S=D  4, the two factors strike a balance so that the average Nusselt number of the upper sphere equals to the lower one. However, when S=D > 4, the effect of mixing gradually exceeds the preheating effect from the lower sphere on the upper sphere, so its average number exceeds that of the lower one, and increases continuously. (3) The local Nusselt number of the lower sphere shows a positive dependence on the spacing and reaches the maximum at the front stagnation point (h ¼ 0 ), which has a similar trend with that of the isothermal single sphere. For h > 135 , the variations of various spacing come to one curve. However, the forward stagnation point of the upper sphere moves to h ¼ 40 due to stronger preheating effect at S=D ¼ 2, and it emerges two confluence points at h  120 and 160 , respectively. In general, the local Nusselt number of both spheres declines gradually as the angle increases. (4) Influenced by the lower sphere, the relative magnitude of the friction and pressure resistance of the upper sphere expresses a reciprocal relationship. The average drag coefficient of the upper sphere increases with the increasing of S=D, but the value of the lower sphere almost does not change. For C DP =C D , the ratio of the lower sphere decreases gradually over the range of 2 6 S=D 6 24, while the ratio of the upper sphere decreases within S=D 6 4 and then keeps rising for 4 < S=D 6 24. (5) The variation of dimensionless temperature with the distance at the centerline shows that when S=D ¼ 24, the lower sphere has little influence on the upper sphere. According to the relationship of Nu  S=D, natural convection heat transfer correlations for accurately predicting the heat transfer rate are proposed in the range of 2 6 S=D 6 24, which can be used in related calculations of engineering applications and academic research. For the lower sphere:

NuL ¼ 12:1417  2:3712  0:9906ðS=DÞ :

For the upper sphere:

NuU ¼ 15:1716  8:8108  0:8746ðS=DÞ : For the whole collection:

NuWhole ¼ 25:5812  9:40  0:8839ðS=DÞ : Conflict of interest The authors declared that there is no conflict of interest. Acknowledgement The National Natural Science Foundation of China (Grant No. 51576189) has supported this work. The National Supercomputing Center (Shenzhen) is also greatly appreciated by the authors, for providing the commercial software, ANSYS Fluent (version 15.0.0), and for their excellent services. References [1] A. Bejan, Convection Heat Transfer, John wiley & sons, 2013. [2] B. Singh, S.K. Dash, Natural convection heat transfer from a finned sphere, Int. J. Heat Mass Transfer 81 (81) (2015) 305–324. [3] O.G. Martynenko, P.P. Khramtsov, Free-Convective Heat Transfer: With Many Photographs of Flows and Heat Exchange, Springer Science & Business Media, 2005. [4] S. Yamoah, E.H. Akaho, N.G. Ayensu, et al., Analysis of fluid flow and heat transfer model for the pebble bed high temperature gas cooled reactor, Res. J. Appl. Sci. Eng. & Tech. 4 (12) (2012) 1659–1666. [5] Y. Lei, W.L. Zhan, New concept for ADS spallation target: gravity-driven dense granular flow target, Sci. China Technol. Sc. 58 (10) (2015) 1705–1711. [6] F. Geoola, A.R.H. Cornish, Numerical simulation of free convective heat transfer from a sphere, Int. J. Heat Mass Transfer 25 (11) (1982) 1677–1687. [7] S.W. Churchill, Comprehensive, theoretically based, correlating equations for free convection from isothermal spheres, Chem. Eng. Commun. 24 (4–6) (1983) 339–352. [8] K. Jafarpur, M.M. Yovanovich, Laminar free convective heat transfer from isothermal spheres: a new analytical method, Int. J. Heat Mass Transfer 35 (9) (1992) 2195–2201. [9] H. Jia, G. Gogos, Laminar natural convection heat transfer from isothermal spheres, Int. J. Heat Mass Transfer 39 (8) (1996) 1603–1615. [10] A. Prhashanna, R.P. Chhabra, Free convection in power-law fluids from a heated sphere, Chem. Eng. Sci. 65 (23) (2010) 6190–6205. [11] K. Kitamura et al., Fluid flow and heat transfer of high-Rayleigh-number natural convection around heated spheres, Int. J. Heat Mass Transfer 86 (2015) 149–157. [12] D.Y. Lee, B.J. Chung, Visualization of natural convection heat transfer on a sphere, Heat Mass Transfer 53 (12) (2017) 3613–3620. [13] M.J. Chamberlain, K.G.T. Hollands, G.D. Raithby, Experiments and theory on natural convection heat transfer from bodies of complex shape, J. Heat Trans. 107 (3) (1985) 624–629. [14] M.M. Yovanovich, On the effect of shape, aspect ratio and orientation upon natural convection from isothermal bodies of complex shape, ASME HTD 82 (1987) 121–129. [15] G.D. Raithby, K.G.T. Hollands, Handbook of Heat Transfer, third ed., McGrawHill Book Company, New York, 1998. [16] A.V. Hassani, K.G.T. Hollands, A simplified method for estimating natural convection heat transfer from bodies of arbitrary shape, ASME National Heat Transfer Conference, Pittsburg, PA, August 8–12, 1987. [17] A.V. Hassani, K.G.T. Hollands, On natural convection heat transfer from threedimensional bodies of arbitrary shape, J. Heat Trans. 111 (2) (1989) 363–371. [18] L. Schouveiler, A. Brydon, T. Leweke, M.C. Thompson, Interactions of the wakes of the two spheres placed side by side, Eur. J. Mech. B-Fluids 23 (2004) 137– 145. [19] E. Pinar, B. Sahin, M. Ozgoren, H. Akilli, Experimental study of flow structures around side-by-side spheres, Ind. Eng. Chem. Res. 52 (40) (2013) 14492– 14503. [20] S. Mukherjee, N. Dasgupta, R.P. Chhabra, Natural Convection from an Array of the two spheres Submerged in Power-Law Fluids, Conference: Chemcon, 2015. [21] S.G. Musong, Z.G. Feng, E.E. Michaelides, S. Mao, Application of a threedimensional immersed boundary method for free convection from single spheres and aggregates, J. Fluid Eng. 138 (4) (2016) 041304. [22] A. Bejan, A.J. Fowler, G. Stanescu, The optimal spacing between horizontal cylinders in a fixed volume cooled by natural convection, Int. J. Heat Mass Transfer 38 (11) (1995) 2047–2055.

J. Zhang et al. / International Journal of Heat and Mass Transfer 134 (2019) 852–865 [23] J. Liu, H. Liu, Q. Zhen, W.Q. Lu, Numerical investigation of the laminar natural convection heat transfer from two horizontally attached horizontal cylinders, Int. J. Heat Mass Transfer 104 (2017) 517–532. [24] J. Liu, H. Liu, Q. Zhen, W.Q. Lu, Laminar natural convection heat transfer from a pair of attached horizontal cylinders set in a vertical array, Appl. Therm. Eng. 115 (2017) 1004–1019. [25] J. Liu, C.J. Zhao, H. Liu, W.Q. Lu, Numerical study of laminar natural convection heat transfer from a hemisphere with adiabatic plane and isothermal hemispherical surface, Int. J. Therm. Sci. 131 (2018) 132–143. [26] J. Zhang, J. Liu, W.Q. Lu, Study on laminar natural convection heat transfer from a hemisphere with uniform heat flux surface, J. Therm. Sci. (2018), https://doi. org/10.1007/s11630-018-1051-y.

865

[27] C. Sasmal, R. Shyam, R.P. Chhabra, Laminar flow of power-law fluids past a hemisphere: momentum and forced convection heat transfer characteristics, Int. J. Heat Mass Transfer 63 (2013) 51–64. [28] P. Mishra, A.K. Tiwari, R.P. Chhabra, Effect of orientation on forced convection heat transfer from a heated cone in Bingham plastic fluids, Int. Commun. Heat Mass 93 (2018) 34–40. [29] R. Shyam, C. Sasmal, R.P. Chhabra, Natural convection heat transfer from two vertically aligned circular cylinders in power-law fluids, Int. J. Heat Mass Transfer 64 (2013) 1127–1152.