- Email: [email protected]
Laminar natural convection in an inclined cylindrical enclosure having finite thickness walls
International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Laminar natural convection in an inclined cylindrical enclosure having finite thickness walls Mikhail A. Sheremet ⇑ Faculty of Mechanics and Mathematics, Tomsk State University, 36, Lenin Avenue, 634050 Tomsk, Russia Institute of Power Engineering, Tomsk Polytechnic University, 30, Lenin Avenue, 634050 Tomsk, Russia
a r t i c l e
i n f o
Article history: Received 3 October 2011 Received in revised form 15 January 2012 Accepted 14 February 2012 Available online 14 April 2012 Keywords: Conjugate heat transfer Natural convection Inclined cylinder Three-dimensional regimes Vector potential Boussinesq approximation
a b s t r a c t Mathematical simulation of unsteady natural convection in an inclined cylinder with heat-conducting walls of finite thickness and a local heat source in conditions of convective heat exchange with an environment has been carried out. Numerical analysis has been based on solution of the convection equations in the dimensionless variables vector potential components, modified vorticity functions, temperature. Particular efforts have been focused on the effects of four types of influential factors such as the Rayleigh number Ra = 104, 5 104, 105, the Prandtl number Pr = 0.7, 7.0, the thermal conductivity ratio k2;1 ¼ 5:7 104 ; 4:3 102 and the inclination angle c = 0, p/6, p/3, p/2 on the velocity and temperature fields. The effect scales of the key parameters on the average Nusselt number have been determined. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Study of natural convection is related to definition of optimum heat transfer modes in various technological systems such as heat pipes [1,2], thermosyphons [3,4], cooling systems of heat-generating components in electronics [5,6], chemical reactors [7]. Correct definition of the most effective conditions of the transport processes evolution in such devices is possible only by multiparameter mathematical simulation of nonstationary convective heat transfer modes [8]. To date, a bundle of experimental and theoretical studies of natural convection regimes in cavities with various shapes [9–26] has been conducted. The majority of the investigations concern the numerical analysis of transport processes in the twodimensional objects both in view of heat-conducting walls effect [11–15], and in case of absence of such influence [16–19]. Thus it is necessary to note essential differences of the obtained results [12,13,15], that is caused by significant thermal lag effect of solid walls. For example, Liaqat and Baytas [13] have numerically analyzed natural convection in an enclosure having both heat-conducting walls of finite thickness and without them. The obtained results reflect strong effect of heat-conducting walls of finite thickness on heat transfer regimes. Sheremet [15] has analyzed the diffusion effects in the conjugate heat and mass transfer problems in ⇑ Address: Faculty of Mechanics and Mathematics, Tomsk State University, 36, Lenin Avenue, 634050 Tomsk, Russia. Tel.: +7 3822 412 462; fax: +7 3822 529 740. E-mail address: [email protected] 0017-9310/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2012.02.046
a wide range of key parameters. It was shown that essential changes of thermohydrodynamic regimes in a cavity and also a significant decrease in the average Nusselt and Sherwood numbers is observed in case of infinitely thin walls. The effect of heat-conducting walls on the velocity and temperature fields has been pointed out for the three-dimensional natural convection regimes in rectangular domains [21–23]. Kuznetsov and Sheremet [21] have shown that for the conjugate Rayleigh–Benard problem a decrease in the thermal conductivity ratio leads both to an increase in temperature in a cavity and to a reduction in the average Nusselt number. Valencia et al. [22] have conducted the experimental and numerical analysis of natural convection in a cubical enclosure with and without heat-conducting walls of finite thickness at 107 6 Ra 6 108 . It has been shown, that in case of the conjugate problem the change in circulation rates and temperature of a fluid in the cavity is observed. Ha and Jung [23] have numerically investigated the effect of the heat-generating cubic conducting body on the flow structure in a vertical cubic enclosure. These authors demonstrated that the presence of the heat-conducting body leads to a significant change in the average Nusselt number. The effect of heat-conducting walls on the velocity and temperature fields has been unfairly neglected at the analysis of threedimensional convective heat transfer in cylindrical enclosures [24–26]. Li et al. [24] have numerically analyzed unsteady threedimensional thermo-hydrodynamic structures in a vertical closed cylinder heated from the side and cooled from above. These authors showed that an increase in temperature difference leads
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3583
Nomenclature Bi = hLr/kq Biot number ffi 1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
Fo ¼ a1 = gbðT hs T 0 ÞL3r Fourier number g acceleration of gravity h heat transfer factor k1 thermal conductivity of solid walls k2 fluid thermal conductivity ki,j = ki/kj thermal conductivity ratio lw solid wall thickness Lr cylinder radius p pressure Pr = m/a2 Prandtl number r radial cylindrical coordinate R dimensionless radial cylindrical coordinate Ra ¼ gbðT hs T 0 ÞL3r =ma2 Rayleigh number t time T temperature T0 initial temperature heat source temperature Ths Vr velocity along the r-axis Vu velocity along the u-axis Vz velocity along the z-axis U dimensionless velocity along the r-axis V dimensionless velocity along the u-axis pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V b ¼ gbðT hs T 0 ÞLr buoyancy velocity W dimensionless velocity along the z-axis z vertical cylindrical coordinate heat source thickness zhs Z vertical cylindrical coordinate dimensionless 2 @2 r2 ¼ 1r @r@ r @r@ þ r12 @@u2 þ @z Laplacian 2
to the formation of unstable heat transfer modes. Specifically, at Ra = 2000 the resulting flow is steady but asymmetric. As Ra reaches 3000 the flow already becomes time periodic and oscillates in a large amplitude. Leong [25] in variables such as the vector potential components and the vorticity vector has numerically solved the three-dimensional Rayleigh–Benard convection equations for a vertical cylinder with infinitely thin walls. He has shown three hydrodynamic patterns in the cylinder depending on the Rayleigh
Fig. 1. A scheme of the system: (1) walls; (2) fluid; and (3) heat source.
2
2
~ 2 ¼ 1 @ R @ þ 12 @ 2 þ @ 2 dimensionless Laplacian r R @R @R @Z R @u Greek symbols a1 thermal diffusivity of solid walls a2 fluid thermal diffusivity ai,j = ai/aj thermal diffusivity ratio b coefficient of volumetric thermal expansion c inclination angle Ds computational time step Dr computational step along radial cylindrical coordinate Du computational step along azimuthal cylindrical coordinate Dz computational step along vertical cylindrical coordinate H dimensionless temperature He environmental dimensionless temperature m kinematic viscosity q2 fluid density s dimensionless time u azimuthal cylindrical coordinate wr, wu, wz vector potential components Wr, Wu, Wz dimensionless vector potential components xr, xu, xz modified vorticity functions Xr, Xu, Xz dimensionless modified vorticity functions Subscripts i, j numbers of the solution domain elements (Fig. 1) avg average e environment hs heat source
number. Cheng et al. [26] investigated influence of the thermal boundary conditions on the three-dimensional convective flow in a vertical cylinder heated from below. These authors revealed that in case of the adiabatic side surface the flow is highly asymmetric and contains multicellular vortices even at symmetric mathematical statement of the problem.
Fig. 2. Variation of the average Nusselt number versus the dimensionless time and the mesh parameters at Ra = 104, Pr = 0.7, k2,1 = 5.7 104, c = p/6.
3584
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 3. Velocity field and isosurface of the vertical velocity component ~ V and temperature field H at k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 for Ra = 104 (a), Ra = 5 104 (b) and Ra = 105 (c).
The objective of the current investigation is mathematical simulation of transient natural convection modes in an inclined cylinder having heat-conducting walls of finite thickness at presence of a local heat source in conditions of convective heat exchange with an environment.
2. Physical and mathematical model The physical model of transient conjugate natural convection in an inclined closed cylinder and the coordinate system are schemat-
ically shown in Fig. 1. The domain of interest consists of a fluid cavity bounded by heat-conducting walls of finite thickness. A heat source located at the bottom of the cavity is kept at constant temperature. The convective heat exchange with an environment is modeled at borders r = Lr and z = z1. The border z = 0 is supposed to be adiabatic. It is assumed in the analysis that the thermophysical properties of the walls material and of the fluid are independent of temperature, and the flow is laminar. The fluid is viscous, heat-conducting, Newtonian, and the Boussinesq approximation is valid. The fluid motion and heat transfer in the cavity are supposed to be
3585
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
three-dimensional. Radiation heat exchange from the heat source and between the walls is neglected in comparison with convection. Taking into account the above-mentioned assumptions, the differential problem of non-stationary conjugate natural convection in an inclined closed cylinder on the basis of conservation laws of mass, momentum and energy can be written as follows [27]: In the fluid cavity
aim of the present investigation is the analysis of a thermal state of the system in conditions of a mutual influence of natural convection inside a vertical cylinder filled with fluid and heat conduction in solid walls. For this purpose we shall enter into consideration satisfying the condition ¼ w i þ w j þ w k, the vector potential w r u z ¼ 0 [29,30]. Then components of the velocity will such as div w be defined as follows:
@ðrV r Þ @V u @ðrV z Þ þ þ ¼ 0; @r @z @u
¼ rot w V 1 @wz @wu @w @w 1 @ðrwu Þ 1 @wr ) Vr ¼ : ; Vu ¼ r z ; Vz ¼ r @u r @r r @u @z @z @r
ð1Þ
2 @V r @V r V u @V r @V r V u þ Vz þ Vr þ @t @r r @u @z r @V 1 @p V 2 u r gbðT T 0 Þ þ m r2 V r 2 2 ¼ r @u q2 @r r
sin c cos u;
ð7Þ The modified vorticity functions will become
ð2Þ
xr ¼ r
@V u @V z ; @z @u
xu ¼
@V z @V r ; @r @z
xz ¼
@V u @V r Vu r : @u @r ð8Þ
@V u @V u V u @V u @V u V r V u þ Vz þ Vr þ þ @t @r r @u @z r V u 2 @V r 1 @p 2 þ gbðT T 0 Þ sin c sin u; þ m r Vu 2 þ 2 ¼ r q2 @ u r @u r ð3Þ @V z @V z V u @V z @V z þ Vz þ Vr þ @t @r r @u @z 1 @p þ mr2 V z þ gbðT T 0 Þ cos c; ¼ q2 @z @T @T V u @T @T þ Vr þ ¼ a2 r2 T; þ Vz @t @r @z r @u
ð4Þ
Z ¼ z=Lr ; s ¼ tV b =Lr ; H ¼ ðT T 0 Þ=ðT hs T 0 Þ; 8 8 8 > < U ¼ V r =V b ; > < Wr ¼ wr =ðV b Lr Þ; > < Xr ¼ xr =V b ; 2 2 2 ~ V ¼ V u =V b ; Wu ¼ wu =ðV b Lr Þ; Xu ¼ xu Lr =V b ; r ¼ Lr r ; > > > : : : W ¼ V z =V b ; Wz ¼ wz =ðV b Lr Þ; Xz ¼ xz =V b :
ð5Þ
The dimensionless equations for the fluid can be written as follows:
R ¼ r=Lr ;
~ 2 Wr Wr 2 R2 r
In the solid walls
@T ¼ a1 r2 T: @t
These functions can be derived by means of the vorticity vector ¼ rot V ~ ¼x ~ uj þ x ~ zk ~ ri þ x like that components x ~ r ; xu ¼ x ~ u ; xz ¼ r x ~ z. xr ¼ rx A mathematical model was formulated in terms of the dimensionless variables the vector potential and the modified vorticity functions by using the non-dimensional variables:
ð6Þ
A transformation of the formulated system of differential equations (1)–(6) to a form eliminating direct search of the pressure field [21,28] is represented as the most reasonable because the
@ Wu ¼ RXr ; @u
~ 2 Wu Wu þ 2 R2 r ~ 2 Wz ¼ Xz ; Rr
@ Wr ¼ R2 Xu ; @u
ð9Þ
ð10Þ ð11Þ
Fig. 4. Isolines of vector potential component Wr and temperature fields H at R = 0.28, k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 for Ra = 104 (a), Ra = 5 104 (b) and Ra = 105 (c).
3586
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 5. Isolines of vector potential component Wu and temperature fields H at u = 0.56p, k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 for Ra = 104 (a), Ra = 5 104 (b) and Ra = 105 (c).
@ Xr @ Xr V @ Xr @ Xr @U @U @U U Xr @V Xu þU þW Xz þ Xr U @R @u @Z @Z @s @R R @ u @Z R rffiffiffiffiffiffi @ X Pr ~ 2 2 @ Xr 2 @H @H u þR r Xr cos c; ¼ sin c sin u Ra R @R R @ u @Z @u ð12Þ
@ Xu @ Xu V @ Xu @ Xu Xr @V Xu @V Xz @V U Xu V Xr þ þ 2 þU þW R @u @s @R @Z R @R R @ u R @Z R R rffiffiffiffiffiffi Xu 2 @ Xr Pr ~ 2 @H @H ¼ r Xu 2 þ 3 þ sin c cos u þ cos c; @Z @R Ra R @u R ð13Þ
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3587
Fig. 6. Isolines of vector potential component Wz and temperature fields H at Z = 0.56, k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 for Ra = 104 (a), Ra = 5 104 (b), Ra = 105 (c).
@ Xz @ Xz V @ Xz @ Xz @W @W @W U Xz @V þU þW Xz þ Xr þU Xu @s @R R @ u @Z R @R @u @Z @R rffiffiffiffiffiffi Pr ~ 2 X 2 @ Xz @H @H r Xz þ 2z sin c cos u R ¼ sin c sin u; @u @R Ra R @R R ð14Þ
@H @H V @H @H 1 ~ 2 H: þU þW þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi r @s @R R @ u @Z Pr Ra
ð15Þ
Energy equation for solid walls:
1 @H ~ 2 ¼ r H: Fo @ s
ð16Þ
Eqs. (9)–(11) have been obtained as a result of substitution of the velocity components (7) in the definition of the modified vorticity functions (8). Eqs. (12)–(14) have been obtained according to the entered modified vorticity functions (8). For example, Eq. (12) has been obtained as follows. Firstly we differentiate Eq. (3) with
3588
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Initial conditions are
Wr ðR; u; Z; 0Þ ¼ 0; Wu ðR; u; Z; 0Þ ¼ 0; Wz ðR; u; Z; 0Þ ¼ 0; Xr ðR; u; Z; 0Þ ¼ 0; Xu ðR; u; Z; 0Þ ¼ 0; Xz ðR; u; Z; 0Þ ¼ 0; H(R, u, Z, 0) = 0 except the temperature for the heat source at which H = 1 during the whole process time. Boundary conditions are: convective heat exchange with an environment is modeled at the walls such as R = 1 and Z = z1/Lr: @@Hn ¼ BiðHe HÞ; at the wall Z = 0 for Eq. (16) the heat insulation condition is set
@H ¼ 0; @Z at the solid–fluid interface R = r1/Lr:
@ Wr ¼ Wu ¼ Wz ¼ 0; @R
(
H1 ¼ H2 ; @ H1 @R
H2 ¼ k2;1 @@R ;
at the solid–fluid interfaces Z = const except the heat source surface: Fig. 7. Variation of the average Nusselt number versus the Rayleigh number and the inclination angle at k2,1 = 5.7 104, s = 100, Pr = 0.7
respect to z. The received equation is multiplied by r and Eq. (4) differentiated with respect to u is subtracted from it. Eqs. (9)–(16) are subjected to the following initial and boundary conditions.
@ Wz Wr ¼ Wu ¼ ¼ 0; @Z
(
H1 ¼ H2 ; @ H1 @Z
H2 ¼ k2;1 @@Z ;
Wz at the heat source surface: Wr ¼ Wu ¼ @@Z ¼ 0; H ¼ 1:0; hs for the heat source ðat lLwr 6 Z 6 lw þz ; 0 6 R 6 rL1r ; 0 6 u < 2pÞ Lr H = 1.0; at u = 0, 2p the periodicity conditions are set, e.g.
Fig. 8. Isolines of vector potential component Wr and temperature fields H at R = 0.28, k2,1 = 5.7 104, Ra = 104, Pr = 0.7, c = p/3 for s = 3 (a), s = 6 (b), s = 20 (c) and s = 70 (d).
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
8 8 8 < Wr ju¼0 ¼ Wr ju¼2p ; < Xr ju¼0 ¼ Xr ju¼2p ; < Hju¼0 ¼ Hju¼2p ; @ H : @@Wur ¼ @@Wur ; : @@Xur ¼ @@Xur ; : @@H ¼ : u @ u u¼0
u¼2p
u¼0
u¼2p
u¼0
u¼2p
The correct boundary conditions at R = 0 for vector potential components and modified vorticity functions and also for temperature have been formulated in accordance with the differential equations. The detailed description of such boundary conditions has been presented by Paasonen [31]. Boundary conditions for the modified vorticity functions are obtained on the basis of the following correlations:
Xr ¼ R
@V @W ; @Z @ u
Xu ¼
@W @U ; @R @Z
Xz ¼
@U @V : V R @u @R
Converting these correlations for the internal boundaries such as R = const, Z = const one can obtain: – Xr ¼ 0; Xu ¼ @W ; Xz ¼ R @V at R ¼ const; @R @R @U – Xr ¼ R @V ; X ¼ ; X ¼ 0 at Z ¼ const: u z @Z @Z
3589
Thus for definition of the modified vorticity functions at internal boundaries it is necessary to determine the velocity field: @W U ¼ 1R @@Wuz @Zu ;
V¼
@ Wr @ Wz ; @Z @R
W¼
1 @ðRWu Þ 1 @ Wr : R @R R @u
Eqs. (9)–(16) with corresponding initial and boundary conditions has been solved by means of finite differences method [21,28,32,33]. Uniform grid in R, u and Z directions was used for all computations. A locally one-dimensional scheme of Samarskii [32] was used to solve Eqs. (12)–(16) numerically. In this scheme a solution of the three-dimensional scheme reduces to a sequential solution of the one-dimensional systems. It is to be noted here that Samarskii monotone approximation was applied to discretization of the convective terms in Eqs. (12)–(15) [33]. An implicit difference scheme was used. An evolutionary term represented a one-sided difference in time and had the first order of accuracy in time step. All derivatives with respect to spatial coordinates were approximated with
Fig. 9. Isolines of vector potential component Wu and temperature fields H at u ¼ 0:56p, k2;1 ¼ 5:7 104 ; Ra ¼ 104 ; Pr ¼ 0:7; c ¼ p=3 for s = 3 (a), s = 6 (b), s = 20 (c) and s = 70 (d).
3590
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 10. Isolines of vector potential component Wz and temperature fields H at Z ¼ 0:56, k2;1 ¼ 5:7 104 ; Ra ¼ 104 ; Pr ¼ 0:7; c ¼ p=3 for s = 3 (a), s = 6 (b), s = 20 (c) and s = 70 (d).
the second order of accuracy in the step along the coordinate. As a result an error of approximation of the partial differential equations and the boundary conditions was O(Ds + (Dr)2 + (D
u)2 + (Dz)2). A computational time step was selected from conditions of a stable algorithm and so Ds = 103 [21,28,34]. It should be noted that analysis of the time step effect was previously
3591
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
8 W W W Wzi;j1;k > U i;j;k ¼ zi;jþ1;k ui;j;kþ12hz ui;j;k1 ; > 2Ri hu > < W W W W V i;j;k ¼ ri;j;kþ12hz ri;j;k1 ziþ1;j;k2hr zi1;j;k ; > > > : W ¼ Wuiþ1;j;k Wui1;j;k þ Wui;j;k Wri;jþ1;k Wri;j1;k : i;j;k 2hr R 2R hu i
Fig. 11. Variation of the average Nusselt number versus the dimensionless time and the Rayleigh number at Pr = 0.7, k2,1 = 5.7 104, c = p/3.
performed with time steps Ds of 0.01, 0.001, 0.0001. It was found that at Ds = 102 increasing oscillations of the average Nusselt number have been observed at high analyzed Rayleigh numbers. The numerical stability has been marked for time steps of 0.001 and 0.0001. Therefore the time step of 0.001 was adopted in the present study. Eqs. (9)–(11) have been discretised by means of seven-points difference scheme, so-called ‘‘cross scheme’’, on the basis of a symmetric approximation of the second derivatives. The obtained system of linear algebraic equations has been solved by the successive over relaxation method. Optimum value of a relaxation parameter has been chosen on the basis of computing experiments. Vorticity boundary conditions have an essential value for numerical simulation of analyzed process. Approximation of these boundary conditions for the internal boundaries such as R = r1/Lr, Z = (lw + zhs)/Lr, Z = (z1 - lw)/Lr with the second order of accuracy on the basis of expansion in a Taylor series can be written as follows
8 Xrnr;j;k ¼ 0; > > < 3W 4W nr1;j;k þW nr2;j;k Xunr;j;k ¼ @W ¼ nr;j;k ; at R ¼ Rnr ¼ r 1 =Lr ; — @R nr;j;k 2hr > > 4V nr1;j;k V nr2;j;k 3V nr;j;k : @V Xznr;j;k ¼ R @R nr;j;k ¼ Rnr 2hr 8 4V i;j;nz þ1 V i;j;nz þ2 3V i;j;nz > 1 1 1 > X ¼ R @V ¼ Ri ; > 2hz @Z i;j;nz1 < ri;j;nz1 lw þ zhs 3U þU 4U i;j;nz1 i;j;nz1 þ2 i;j;nz1 þ1 — @U X ¼ @Z i;j;nz ¼ ; at Z ¼ Z nz1 ¼ Lr ; > 2hz 1 > ui;j;nz1 > : Xzi;j;nz1 ¼ 0 8 @V 3V i;j;nz 4V i;j;nz 1 þV i;j;nz 2 > 2 2 2 > ; > Xri;j;nz2 ¼ R @Z i;j;nz2 ¼ Ri 2hz < z 1 lw 4U U 3U i;j;nz 1 i;j;nz 2 i;j;nz — 2 2 2 Xui;j;nz2 ¼ @U ¼ ; atZ ¼ Z nz2 ¼ Lr ; > 2h @Z i;j;nz z 2 > > : Xzi;j;nz2 ¼ 0 where i, j, k are mesh point coordinates. For definition of the velocity field the second-order central differences relations are used for spatial derivatives such as
i
The concise description of numerical simulation of boundary conditions at a cylinder axis in the present work is impossible because the numerical approximation becomes complicated not only owing to the analysis of a nonsymmetric case, but also owing to the selected variables such as vector potential components and modified vorticity functions. The detailed description of the chosen method has been completely presented by Paasonen [31]. It should be noted that for realization of the coordinate-splitting method (a locally one-dimensional scheme of Samarskii [32]) the correct sequence in a choice of coordinates is important. Firstly, the solution along circles with the periodicity conditions at u = 0 and u = 2p is carried out, then the solution along all internal grid lines parallel to an axis of the cylinder, except for this axis, with conditions at the bases of cylinder is conducted. Further on at the cylinder axis the equation obtained by means of the coordinate-splitting method is solved. At the third stage there is a solution of the system of the one-dimensional equations along the radiuses of the cylinder by means of the method presented in [31]. The method of solution has been tested on different meshes. Fig. 2 shows sensitivity of the average Nusselt number at the heat R 0:88 R 2p @ H 1 source surface Nuavg ¼ 1:76 j @Z jZ¼0:24 du dR to the mesh 0 p 0 parameters with increase in time at Ra = 104, Pr = 0.7, c = p/6, k2,1 = 5.7 104. On the ground of optimization of calculation accuracy for 0 6 s 6 100 the computational mesh such as 50 50 50 has been selected. There were 10 nodes along the R-coordinate inside the solid walls. 3. Results and discussion Numerical investigation of the boundary value problem (9)– (16) has been carried out at following values of dimensionless complexes: Ra = 104, 5 104, 105; Pr = 0.7, 7.0; 0 6 s 6 100; c = 0, p/6, p/3, p/2; He = 1; k2,1 = 5.7 104 – Bi = 4.35 102, Fo = 3.5 104; and k2,1 = 4.3 102 – Bi = 3.33, Fo = 9.38 106. Analysis has been carried out for the following geometrical parameters: z1/Lr = 2, lw/Lr = 0.12, r1/Lr = 0.88, zhs/Lr = 0.12. Particular efforts have been focused on the effects of four types of influential factors such as the Rayleigh and Prandtl numbers, the thermal conductivity ratio and an inclination angle on the thermohydrodynamic regimes. 3.1. Effect of Rayleigh number Velocity and temperature fields and isosurfaces of the vertical velocity component at k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 and different values of the Rayleigh number are presented in Fig. 3. This figure reflects formation of the vortex structures in the gas cavity having the inclination angle c = p/3. The motion paths of gas particles crossing the cylinder axis demonstrate the presence of inclination angle. More detailed analysis of distributions of thermohydrodynamic parameters will be presented for specific sections of the domain of interest. Figs. 4–6 show the isolines of the vector potential components and the isotherms for various cross-sections at k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 and different values of the Rayleigh number Ra = 104, 5 104, 105. The increase in the Rayleigh number for the considered range is reflected in an insignificant modification of the vector potential component Wr isolines (Fig. 4). Two convective cells are formed
3592
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 12. Velocity field and isosurfaces of the vertical velocity component ~ V and temperature field H at Ra = 104, k2,1 = 5.7 104, s = 100, c = p/6 for Pr = 0.7 (a) and Pr = 7.0 (b).
in the cross-section R = 0.28. These cells differ in the direction of fluid motion and reflect formation of a thermal plume symmetrically with respect to axis u = p. It is necessary to note the displacement of a convective cells core to the axis u = p with the increase in the Rayleigh number that is caused by an intensification of the transport processes of momentum and energy in this zone. Appearance of a conductive heat transfer zone near the heat source surface is obvious. It is confirmed by the isotherms position. The increase in the Rayleigh number also leads to a decrease in the qffiffiffiffi
Pr thermal boundary layer thickness O 4 Ra [35] above the heat source. In all probability, the formation of a thermal plume and convective cells symmetrically with respect to the axis u = p, is caused by the presence of the cylinder inclination angle such as c = p/3. Interaction of the low temperature front spreading from the border Z = z1/Lr deep into the domain of interest, with the
thermal plume is observed also. Thus an increase in the temperature differences, owing to a rise in a heat source temperature, leads to a decrease in the sizes of the zone that is subjected to the influence of the low temperature. The increase in the Rayleigh number results not only in an increment of the temperature for the upper layers of the analyzed object, but also in an increase in a thermal plume thickness that is reflected in formation of a significant cross-section energy flow. Analyzing changes of the vector potential component Wu for the cross-section u = 0.56p (Fig. 5) it is possible to note an increase in sizes and intensity of a circulating flow located in the angular zone of the domain of interest. The latter leads to distortion of the basic convective cell and to reduction of the fluid circulation velocities in it. Displacement of a core of the basic vortex is caused by both a development of a recirculation flow described above and a formation of a thermal plume over the heat source close to the solid walls.
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3593
Fig. 13. Isolines of vector potential component Wr and temperature fields H at R = 0.28, k2,1 = 5.7 104, Ra = 104, c = p/6, s = 100 for Pr = 0.7 (a) and Pr = 7.0 (b).
Distribution of the isotherms in cross-section u = 0.56p as well as at R = 0.28 reflects both an increase in the temperature for the upper layers of the domain of interest that interferes with penetration of the low temperature front and a drift of convective cell core to zone 0.5 < R < 0.7. The presence of thermal spots in the cavity characterizes an essential three-dimensionality of the analyzed process. Distributions of the vector potential component Wz and the isotherms for the cross-section Z = 0.56 (Fig. 6) confirm the presence of symmetry with respect to the axis u = p observable in Fig. 4, and also reflect both an intensification of convective flow at the increase in the Rayleigh number and a heating of the fluid cavity. The analysis of the influence of the Rayleigh number and the inclination angle on the average Nusselt number at the heat source surface (Fig. 7) has been conducted. An increase in the generalized heat transfer coefficient at the increment in the Rayleigh number is observed. The increase in the inclination angle leads to convective heat transfer intensification owing to an increase in the temperature gradient at the heat source surface. The nonmonotonic changes in the average Nusselt number at a variation of the inclination angle a reflects the presence of an optimum position of the cylinder with the maximum value of Nuavg. In the analyzed range of the Rayleigh number it is possible to indicate two intervals with the maximum magnitude of the generalized heat transfer coefficient at various values of the inclination angle of the cylinder such as at 104 6 Ra 6 3 104 we have c = p/3 and at 3 104 < Ra 6 105 we have c = p/6. The presented average Nusselt number curves in Fig. 7 have been defined as correlations that are very convenient in carrying out applied calculations:
Nuavg ¼ 0:188 Ra0:315 at c ¼ 0; Nuavg ¼ 0:242 Ra0:305 at c ¼ p=6; Nuavg ¼ 0:323 Ra0:277 at c ¼ p=3; Nuavg ¼ 0:34 Ra0:265 at c ¼ p=2:
3.2. Effect of the transient factor
Fig. 14. Isolines of vector potential component Wu and temperature fields H at u = 0.56p, k2,1 = 5.7 104, Ra = 104, c = p/6, s = 100 for Pr = 0.7 (a) and Pr = 7.0 (b).
The transient factor in conjugate heat transfer problems plays an essential role [21,28] as it defines not only dynamics of the velocity and temperature fields owing to formation, evolution and dissipation of the vortex structures in a fluid cavity, but also a thermal lag of solid walls in conditions of heat exchange with an environment.
3594
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 15. Isolines of vector potential component Wz and temperature fields H at Z = 0.56, k2,1 = 5.7 104, Ra = 104, c = p/6, s = 100 for Pr = 0.7 (a) and Pr = 7.0 (b).
Fig. 8 demonstrates isolines of the vector potential component
Wr and temperature fields H in the cross-section R = 0.28 at k2,1 = 5.7 104, Ra = 104, Pr = 0.7, c = p/3 for different values of dimensionless time. Two convective cells of small intensity are formed symmetrically with respect to axis u = p at the initial time s = 3 in cross-section R = 0.28 that is caused by conductive heat transfer in the cavity (isotherms are parallel to the heat source surface). The increase in dimensionless time leads to formation and evolution of a thermal plume which intensifies convective heat transfer in the cavity. The sizes of convective cells are increased while cores of the vortexes are displaced to a center layer of the cylinder. Penetration of a low temperature front from the border Z = z1/Lr deep into the cavity is observed jointly with evolution of the thermal plume. In the cross-section u = 0.56p (Fig. 9) the dynamics of formation and evolution of the hydrodynamic structures also reflects the basic moments described above. The increase in the dimensionless time s leads to an intensification of the basic convective cell formed directly above the heat source. An appearance of a vortex is observed in the zone 1 < Z < 1.8 owing to cooling of the top wall of the domain of interest. This vortex loses an internal energy and decreases in sizes at the increase in dimensionless time. Distribution of isotherms in the fluid cavity at s P 20 is defined not only by influence of the buoyancy force, but also by dynamics of con-
Fig. 16. Variation of the average Nusselt number versus the dimensionless time and the Prandtl number at Ra = 104, k2,1 = 5.7 104, c = p/6.
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
ductive heat transfer in the solid walls of the cylinder. Presence of the inclination angle such as c = p/3 and heating of the solid walls lead to formation of the thermal plume near to the side surface of the cylinder. The modification of the isolines of the vector potential component Wz in the cross-section Z = 0.56 (Fig. 10) reflects an essential nonstationarity of the analyzed process. Formation of two groups of low-intensity convective cells differing by sizes is observed at the initial time. Disintegration of the large vortexes to the finer vortexes and reorientation of formed hydrodynamic structures are observed at the increase in the dimensionless time. There is an increase in sizes of the central vortexes group at s = 20; and by the moment s = 70 this two-vortical system becomes dominating. It is necessary to note, that the hydrodynamic
3595
structures formed at the moments s = 3 and s = 70 are similar to each other. A basic difference consists in orientation and intensity of these vortexes. The temperature field varies insignificantly. An increase in the thickness of thermal plume located near to the side surface of the cylinder in the zone u = p is observed. Dynamics of the generalized heat transfer coefficient at the heat source surface at various values of the Rayleigh number is presented in Fig. 11. Change of the average Nusselt number can be broken into two levels. 0 < s < 20 is an initial level reflecting presence of oscillations in a function of Nuavg(s). 20 < s < 100 is a level of a thermal stabilization of natural convection. It is necessary to note, that at large Rayleigh numbers ðRa P 105 Þ stabilization is delayed owing to essential interaction of the temperature fields of solid walls and the fluid cavity.
Fig. 17. Velocity field and isosurfaces of the vertical velocity component ~ V and temperature field H at Ra = 105, Pr = 0.7, c = p/6, s = 100 for k2,1 = 5.7 104 (a) and k2,1 = 4.3 102 (b).
3596
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3.3. Effect of the Prandtl number Further the effect of fluid heat-transfer properties on the flow and heat transfer modes will be analyzed. The increase in the Prandtl number in 10 times leads to attenuation of convective heat transfer in the cavity (Figs. 12–15,b), that is caused by an increase in a kinematic viscosity coefficient and an internal friction force. The latter is reflected as well in both velocity field (Fig. 12) and configurations of isolines of the vector potential components e.g. Wz (Fig. 15). Slight displacement of the convective cells cores at Pr = 7.0 in the cross-section R = 0.28 (Fig. 13) is caused by more intensive diffusion of perturbations in the cavity from solid walls. An increase in sizes of a recirculation zone located in an angular area of the fluid cavity is observed in the cross-section u = 0.56p (Fig. 14). The latter is related also to heat conduction dynamics in solid walls. It is possible to note more intensive heating of the domain of interest at Pr = 7.0 (Fig. 12) that is caused by both hydrodynamics, namely, small flow velocities, presence of recirculation zones reflecting formation of high temperature areas, and thermodynamics of the conjugate problems, namely, interaction of conductive heat transfer in solid walls and natural convection in the fluid cavity. The latter leads to a decrease in an intensity of the low temperature penetration from the border Z = z1/Lr deep into the domain of interest (Figs. 12 and 14), and also to an increase in a thermal plume thickness (Figs. 12–15). It is necessary to note formation of the thermal spots (Figs. 13 and 14), that is caused by the three-dimensionality of the analyzed process. Change dynamics of the generalized heat transfer coefficient versus the Prandtl number is presented in Fig. 16. An increase in the fluid viscosity is reflected in a delaying of the thermal stabilization zone. At initial time s < 47 a maximum intensity of convective heat transfer is observed at Pr = 7.0, and at s > 47 vice versa. 3.4. Effect of the thermal conductivity ratio Presence of heat-conducting walls of finite thickness can essentially modify the temperature field in the fluid cavity and flow regimes [12,15]. The solid walls can play a role as an original buffer zone owing to a thermal lag of the solid material. For a concrete time interval walls can accumulate energy, reflecting indirect influence of an environment on the processes proceeding in the cavity. Figs. 17–20 show the velocity and temperature fields and also the isolines of the vector potential components in various crosssections, corresponding to convective heat transfer mode
Fig. 19. Isolines of vector potential component Wu and temperature fields H at u = 0.56p, Ra = 105, Pr = 0.7, c = p/6, s = 100 for k2,1 = 5.7 104 (a) and k2,1 = 4.3 102 (b).
Fig. 18. Isolines of vector potential component Wr and temperature fields H at R = 0.28, Ra = 105, Pr = 0.7, c = p/6, s = 100 for k2,1 = 5.7 104 (a) and k2,1 = 4.3 102 (b).
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3597
Fig. 20. Isolines of vector potential component Wz and temperature fields H at Z = 0.56, Ra = 105, Pr = 0.7, c = p/6, s = 100 for k2,1 = 5.7 104 (a) and k2,1 = 4.3 102 (b).
Ra = 105, Pr = 0.7, s = 100, c = p/6, at different values of the thermal conductivity ratio k2,1 = 5.7 104, 4.3 102. The increase in the thermal conductivity ratio is caused with a decrease in the thermal conductivity of the solid material. Transition from k2,1 = 5.7 104 to k2,1 = 4.3 102 leads to essential modification of the vector potential component Wz in the cross-section Z = 0.56 (Fig. 20). Combination of the vortical structures with formation of the global circulation defining mass transfer, momentum transport and energy transport in the cavity is observed. The flow intensity is increased, though slight decrease in the fluid velocities is observed in a vertical direction (Fig. 18). The most essential changes concern the temperature fields. The temperature in the fluid cavity is decreased at k2,1 = 4.3 102 (Fig. 17) that is related to a diminution of the thermal conductivity of solid walls and, accordingly, with low speed of the energy transport in solid walls both from the heat source and from an environment. The latter is evidently presented in Fig. 19 where the isotherm corresponding to the dimensionless temperature H = 0.2 at k2,1 = 4.3 102 is inside of the top wall and in case of k2,1 = 5.7 104 it is at the solid–fluid interface. It is necessary to note loss of stability of the thermal plume above the energy source (Fig. 19). It reflects formation of transition thermohydrodynamic regimes. More intensive cooling of the upper wall in case of a heat-conductive material leads to an increase in a temperature gradient at the internal solid–fluid interface. The latter leads to an increase in convective flow intensity in the cavity.
Fig. 21. Variation of the average Nusselt number versus the Rayleigh number and the thermal conductivity ratio at Pr = 0.7, c = p/6, s = 100.
3598
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 22. Isolines of vector potential component Wr and temperature fields H at R = 0.28, Ra = 104, Pr = 0.7, k2,1 = 5.7 104, s = 100 for c = 0 (a) and c = p/2 (b).
Symmetry of the vector potential components isolines and isotherms with respect to an axis u = p is conserved (Fig. 18), that is caused by both the inclination angle of the cylinder and with indirect influence of an environment. Variation of the average Nusselt number at the heat source surface depending on the Rayleigh number in a range 104 6 Ra 6 105 and the thermal conductivity ratio is presented in Fig. 21. The increase in k2,1 leads to formation of two zones reflecting various behavior of the generalized heat transfer coefficient. The increase in the thermal conductivity ratio leads to an increase in the average Nusselt number at 104 6 Ra < 5:6 104 owing to rising of the temperature gradient that is caused by the growth of energy removal in solid walls. The opposite situation caused by an increase in the role of the hydrodynamic component of the analyzed process is observed in case of 5:6 104 < Ra 6 105 . An increase in the buoyancy force leads to an intensive vertical energy transport in the cavity i.e. the near-surface layers close to the heat source are considerably heated. 3.5. Effect of the inclination angle
Fig. 23. Isolines of vector potential component Wu and temperature fields H at u = 0.56p, Ra = 104, Pr = 0.7, k2,1 = 5.7 104, s = 100 for c = 0 (a) and c = p/2 (b).
Figs. 22–24 present the isolines of vector potential components and the temperature fields in the various cross-sections, corresponding to the convective heat transfer regime at Ra = 104, Pr = 0.7, k2,1 = 5.7 104, s = 100 for various values of the inclination angle c = 0 [34] and c = p/2 . Presence of the inclination angle such as c = p/2 leads to displacement of a thermal plume to an axis u = p and in an intensification of convective flow in the cavity (Fig. 22). Circulating flows reflecting a fluid motion in an opposite direction in comparison with the regime at c = 0 [34] are shaped in the cross-section u = 0.56p (Fig. 23). Clear differences in hydrodynamic structures are displayed in the cross-section Z = 0.56 (Fig. 24) reflecting a system rotation by the angle c = p/2. The temperature fields also differ. There is more intensive cooling of the cavity from an environment in the crosssection u = 0.56p at c = 0. Variations of the average Nusselt number versus the dimensionless time at different values of the inclination angle are presented in Fig. 25. It is necessary to note on the one hand an essential increase in Nuavg at presence of an inclination angle of the cylinder and on the other hand a smoother stabilization of the process. While reorganization of the thermodynamic structures at the moment s = 58 at c = 0 defining presence of a breakpoint in dependence Nuavg(s) is observed. The increase in an inclination
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3599
Fig. 24. Isolines of vector potential component Wz and temperature fields H at Z = 0.56, Ra = 104, Pr = 0.7, k2,1 = 5.7 104, s = 100 for c = 0 (a), c = p/2 (b).
angle leads to an increase in the generalized heat transfer coefficient in a zone s < 10, subsequently this dynamics is broken, and also is reflected in faster stabilization of the analyzed process. 4. Conclusions
Fig. 25. Variation of the average Nusselt number versus the dimensionless time and the inclination angle of the cylinder at Ra = 104, Pr = 0.7, k2,1 = 5.7 104.
Transient three-dimensional natural convection in an inclined closed tube with heat-conducting walls of finite thickness at presence of the internal heat source of the constant temperature in conditions of convective heat exchange with an environment has been numerically solved. Distributions of isolines of the vector potential components and the temperature fields describing features of the analyzed flow and heat transfer regimes at Ra = 104, 5 104, 105, Pr = 0.7, 7.0, 0 6 s 6 100, c = 0, p/6, p/3, p/2, k2,1 = 5.7 104, 4.3 102 have been obtained. A detailed analysis of influence of the Rayleigh number and the Prandtl number, of the dimensionless time, of the inclination angle of the tube and of the thermal conductivity ratio both on thermohydrodynamic fields and on the average Nusselt number at the heat source surface has been conducted. It has been ascertained that in the analyzed range of the Rayleigh number it is possible to indicate two intervals with a maximum magnitude of the generalized heat transfer coefficient at various values of the inclination angle of the tube such as at 104 6 Ra 6 3 104 we have c = p/3 and at 3 104 < Ra 6 105 we have c = p/6 (Fig. 7). It has been shown that variation of the average Nusselt number can be broken into two levels (Fig. 11).
3600
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
0 < s < 20 is an initial level reflecting presence of oscillations in a function Nuavg(s). 20 < s < 100 is a level of thermal stabilization of natural convection. The increase in the Prandtl number is reflected in a delaying of the thermal stabilization zone (Fig. 16). At initial time s < 47 maximum intensity of convective heat transfer is observed at Pr = 7.0, and at s > 47 maximum intensity of convective heat transfer is observed at Pr = 0.7. The various mechanisms defining change of the average Nusselt number versus the thermal conductivity ratio have been ascertained (Fig. 21). The increase in k2,1 leads to an increase in Nuavg at 104 6 Ra < 5:6 104 owing to a prevailing influence of thermal component of the analyzed process. An opposite situation caused by an increase in a role of the hydrodynamic component of natural convection is observed in case of 5:6 104 < Ra 6 105 . An essential increase in Nuavg and a smoother stabilization of the process at presence of the inclination angle of the cylinder have been shown (Fig. 25). Acknowledgements This work was performed according to the Russian Federal Purpose Program during 2009–2013 years (state contract No. P 357) and was supported by the Grants Council (under the President of the Russian Federation), Grant No. MR-5652.2012.8. References [1] G. Carbajal, C.B. Sobhan, G.P. ‘‘Bud’’ Peterson, D.T. Queheillalt, Haydn N.G. Wadley, A quasi-3D analysis of the thermal performance of a flat heat pipe, Int. J. Heat Mass Transfer 50 (2007) 4286–4296. [2] L.L. Vasiliev, Heat pipes in modern heat exchangers, Appl. Therm. Eng. 25 (2005) 1–19. [3] G.V. Kuznetsov, M.A. Al-Ani, M.A. Sheremet, Numerical analysis of convective heat transfer in a closed two-phase thermosyphon, J. Eng. Thermophys. 20 (2011) 201–210. [4] Gilles Desrayaud, Alberto Fichera, Manuel Marcoux, Numerical investigation of natural circulation in a 2D-annular closed-loop thermosyphon, Int. J. Heat Fluid Flow 27 (2006) 154–166. [5] L.A. Florio, A. Harnoy, Combination technique for improving natural convection cooling in electronics, Int. J. Therm. Sci. 46 (2007) 76–92. [6] H. Oprins, J. Danneels, B. Van Ham, B. Vandevelde, M. Baelmans, Convection heat transfer in electrostatic actuated liquid droplets for electronics cooling, Microelectron. J. 39 (2008) 966–974. [7] A.N. Campbell, S.S.S. Cardoso, A.N. Hayhurst, A comparison of measured temperatures with those calculated numerically and analytically for an exothermic chemical reaction inside a spherical batch reactor with natural convection, Chem. Eng. Sci. 62 (2007) 3068–3082. [8] Y. Jaluria, Design and Optimization of Thermal Systems, McGraw-Hill, New York, 1998. p. 626. [9] V.I. Polezhaev, Free convection in the internal problem: results and prospects, J. Eng. Phys. Thermophys. 69 (1996) 676–687. [10] M.V. Diomidov, M.I. Nizovtsev, V.I. Terekhov, Ventilation of window interplane cavity aimed at a higher temperature of the inner pane, Therm. Sci. 6 (2002) 15–22. [11] D.A. Kaminski, C. Prakash, Conjugate natural convection in a square enclosure effect of conduction on one of the vertical walls, Int. J. Heat Mass Transfer 29 (1986) 1979–1988.
[12] A. Liaqat, A.C. Baytas, Numerical comparison of conjugate and non-conjugate natural convection for internally heated semi-circular pools, Int. J. Heat Fluid Flow 22 (2001) 650–656. [13] A. Liaqat, A.C. Baytas, Conjugate natural convection in a square enclosure containing volumetric sources, Int. J. Heat Mass Transfer 44 (2001) 3273– 3280. [14] G.V. Kuznetsov, M.A. Sheremet, Conjugate heat transfer in an enclosure under the condition of internal mass transfer and in the presence of the local heat source, Int. J. Heat Mass Transfer 52 (2009) 1–8. [15] M.A. Sheremet, The influence of cross effects on the characteristics of heat and mass transfer in the conditions of conjugate natural convection, J. Eng. Thermophys. 19 (2010) 119–127. [16] Yasin Varol, Hakan F. Oztop, Ahmet Koca, Filiz Ozgen, Natural convection and fluid flow in inclined enclosure with a corner heater, Appl. Therm. Eng. 29 (2009) 340–350. [17] Chengwang Lei, S.W. Armfield, J.C. Patterson, Unsteady natural convection in a water-filled isosceles triangular enclosure heated from below, Int. J. Heat Mass Transfer 51 (2008) 2637–2650. [18] Fu-Yun Zhao, Di Liu, Guang-Fa Tang, Natural convection in an enclosure with localized heating and salting from below, Int. J. Heat Mass Transfer 51 (2008) 2889–2904. [19] V.I. Terekhov, V.V. Terekhov, Heat transfer in a high vertical enclosure with fins attached to one of the side walls, High Temp. 44 (2006) 436–441. [20] L.A. Moiseeva, S.G. Cherkasov, Steady-state natural-convection heat transfer in a cylindrical tank under conditions of uniform heat supply and simultaneous heat removal through local sinks, High Temp. 35 (1997) 553–558. [21] G.V. Kuznetsov, M.A. Sheremet, The Rayleigh–Benard instability in an enclosure having finite thickness walls, J. Phys: Conf. Ser. 216 (2010) 1–15. [22] L. Valencia, J. Pallares, I. Cuesta, F. Xavier Grau, Turbulent Rayleigh–Benard convection of water in cubical cavities: a numerical and experimental study, Int. J. Heat Mass Transfer 50 (2007) 3203–3215. [23] M.Y. Ha, M.J. Jung, A numerical study on three-dimensional conjugate heat transfer of natural convection and conduction in a differentially heated cubic enclosure with a heat-generating cubic conducting body, Int. J. Heat Mass Transfer 43 (2000) 4229–4248. [24] Y.H. Li, K.C. Lin, T.F. Lin, Computation of unstable liquid metal convection in a vertical closed cylinder heated from the side and cooled from above, Numer. Heat Transfer Part A: Appl. 32 (1997) 289–309. [25] S.S. Leong, Numerical study of Rayleigh–Benard convection in a cylinder, Numer. Heat Transfer Part A: Appl. 41 (2002) 673–683. [26] T.C. Cheng, Y.H. Li, T.F. Lin, Effects of thermal boundary condition on buoyancy driven transitional air flow in a vertical cylinder heated from below, Numer. Heat Transfer Part A: Appl. 37 (2000) 917–936. [27] L.G. Loicyanski, Classical Fluid Mechanics, Drofa, Moscow, 2003 (p. 840). [28] G.V. Kuznetsov, M.A. Sheremet, Conjugate natural convection in an enclosure with local heat sources, Comput. Therm. Sci. 1 (2009) 341–360. [29] K. Aziz, J.D. Hellums, Numerical solution of three-dimensional equations of motion for laminar natural convection, Phys. Fluids 10 (1967) 314–324. [30] G.J. Hirasaki, J.D. Hellums, A general formulation of the boundary conditions on the vector potential in three dimensional hydrodynamics, Quart. Appl. Math. 16 (1968). [31] V.I. Paasonen, Boundary conditions of high-order accuracy at the poles of curvilinear coordinate systems, Russ. J. Numer. Anal. Math. Model. 14 (1999) 369–382. [32] A.A. Samarskii, Theory of Difference Schemes, Nauka, Moscow, 1977. 656 p. [33] V.M. Paskonov, V.I. Polezhaev, L.A. Chudov, Numerical Modelling of Heat and Mass Exchange Processes, Nauka, Moscow, 1984. 288 p. [34] M.A. Sheremet, Three-dimensional conjugate natural convection in a vertical cylinder under heat transfer to the surroundings, Fluid Dyn. 46 (2011) 647– 657. [35] B. Gebhart, Y. Jaluria, R.L. Mahajan, B. Sammakia, Buoyancy – Induced Flows and Transport, Hemisphere Publishing Corporation, New York, 1988. 1001 p.
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Laminar natural convection in an inclined cylindrical enclosure having finite thickness walls Mikhail A. Sheremet ⇑ Faculty of Mechanics and Mathematics, Tomsk State University, 36, Lenin Avenue, 634050 Tomsk, Russia Institute of Power Engineering, Tomsk Polytechnic University, 30, Lenin Avenue, 634050 Tomsk, Russia
a r t i c l e
i n f o
Article history: Received 3 October 2011 Received in revised form 15 January 2012 Accepted 14 February 2012 Available online 14 April 2012 Keywords: Conjugate heat transfer Natural convection Inclined cylinder Three-dimensional regimes Vector potential Boussinesq approximation
a b s t r a c t Mathematical simulation of unsteady natural convection in an inclined cylinder with heat-conducting walls of finite thickness and a local heat source in conditions of convective heat exchange with an environment has been carried out. Numerical analysis has been based on solution of the convection equations in the dimensionless variables vector potential components, modified vorticity functions, temperature. Particular efforts have been focused on the effects of four types of influential factors such as the Rayleigh number Ra = 104, 5 104, 105, the Prandtl number Pr = 0.7, 7.0, the thermal conductivity ratio k2;1 ¼ 5:7 104 ; 4:3 102 and the inclination angle c = 0, p/6, p/3, p/2 on the velocity and temperature fields. The effect scales of the key parameters on the average Nusselt number have been determined. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction Study of natural convection is related to definition of optimum heat transfer modes in various technological systems such as heat pipes [1,2], thermosyphons [3,4], cooling systems of heat-generating components in electronics [5,6], chemical reactors [7]. Correct definition of the most effective conditions of the transport processes evolution in such devices is possible only by multiparameter mathematical simulation of nonstationary convective heat transfer modes [8]. To date, a bundle of experimental and theoretical studies of natural convection regimes in cavities with various shapes [9–26] has been conducted. The majority of the investigations concern the numerical analysis of transport processes in the twodimensional objects both in view of heat-conducting walls effect [11–15], and in case of absence of such influence [16–19]. Thus it is necessary to note essential differences of the obtained results [12,13,15], that is caused by significant thermal lag effect of solid walls. For example, Liaqat and Baytas [13] have numerically analyzed natural convection in an enclosure having both heat-conducting walls of finite thickness and without them. The obtained results reflect strong effect of heat-conducting walls of finite thickness on heat transfer regimes. Sheremet [15] has analyzed the diffusion effects in the conjugate heat and mass transfer problems in ⇑ Address: Faculty of Mechanics and Mathematics, Tomsk State University, 36, Lenin Avenue, 634050 Tomsk, Russia. Tel.: +7 3822 412 462; fax: +7 3822 529 740. E-mail address: [email protected] 0017-9310/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2012.02.046
a wide range of key parameters. It was shown that essential changes of thermohydrodynamic regimes in a cavity and also a significant decrease in the average Nusselt and Sherwood numbers is observed in case of infinitely thin walls. The effect of heat-conducting walls on the velocity and temperature fields has been pointed out for the three-dimensional natural convection regimes in rectangular domains [21–23]. Kuznetsov and Sheremet [21] have shown that for the conjugate Rayleigh–Benard problem a decrease in the thermal conductivity ratio leads both to an increase in temperature in a cavity and to a reduction in the average Nusselt number. Valencia et al. [22] have conducted the experimental and numerical analysis of natural convection in a cubical enclosure with and without heat-conducting walls of finite thickness at 107 6 Ra 6 108 . It has been shown, that in case of the conjugate problem the change in circulation rates and temperature of a fluid in the cavity is observed. Ha and Jung [23] have numerically investigated the effect of the heat-generating cubic conducting body on the flow structure in a vertical cubic enclosure. These authors demonstrated that the presence of the heat-conducting body leads to a significant change in the average Nusselt number. The effect of heat-conducting walls on the velocity and temperature fields has been unfairly neglected at the analysis of threedimensional convective heat transfer in cylindrical enclosures [24–26]. Li et al. [24] have numerically analyzed unsteady threedimensional thermo-hydrodynamic structures in a vertical closed cylinder heated from the side and cooled from above. These authors showed that an increase in temperature difference leads
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3583
Nomenclature Bi = hLr/kq Biot number ffi 1 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
Fo ¼ a1 = gbðT hs T 0 ÞL3r Fourier number g acceleration of gravity h heat transfer factor k1 thermal conductivity of solid walls k2 fluid thermal conductivity ki,j = ki/kj thermal conductivity ratio lw solid wall thickness Lr cylinder radius p pressure Pr = m/a2 Prandtl number r radial cylindrical coordinate R dimensionless radial cylindrical coordinate Ra ¼ gbðT hs T 0 ÞL3r =ma2 Rayleigh number t time T temperature T0 initial temperature heat source temperature Ths Vr velocity along the r-axis Vu velocity along the u-axis Vz velocity along the z-axis U dimensionless velocity along the r-axis V dimensionless velocity along the u-axis pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V b ¼ gbðT hs T 0 ÞLr buoyancy velocity W dimensionless velocity along the z-axis z vertical cylindrical coordinate heat source thickness zhs Z vertical cylindrical coordinate dimensionless 2 @2 r2 ¼ 1r @r@ r @r@ þ r12 @@u2 þ @z Laplacian 2
to the formation of unstable heat transfer modes. Specifically, at Ra = 2000 the resulting flow is steady but asymmetric. As Ra reaches 3000 the flow already becomes time periodic and oscillates in a large amplitude. Leong [25] in variables such as the vector potential components and the vorticity vector has numerically solved the three-dimensional Rayleigh–Benard convection equations for a vertical cylinder with infinitely thin walls. He has shown three hydrodynamic patterns in the cylinder depending on the Rayleigh
Fig. 1. A scheme of the system: (1) walls; (2) fluid; and (3) heat source.
2
2
~ 2 ¼ 1 @ R @ þ 12 @ 2 þ @ 2 dimensionless Laplacian r R @R @R @Z R @u Greek symbols a1 thermal diffusivity of solid walls a2 fluid thermal diffusivity ai,j = ai/aj thermal diffusivity ratio b coefficient of volumetric thermal expansion c inclination angle Ds computational time step Dr computational step along radial cylindrical coordinate Du computational step along azimuthal cylindrical coordinate Dz computational step along vertical cylindrical coordinate H dimensionless temperature He environmental dimensionless temperature m kinematic viscosity q2 fluid density s dimensionless time u azimuthal cylindrical coordinate wr, wu, wz vector potential components Wr, Wu, Wz dimensionless vector potential components xr, xu, xz modified vorticity functions Xr, Xu, Xz dimensionless modified vorticity functions Subscripts i, j numbers of the solution domain elements (Fig. 1) avg average e environment hs heat source
number. Cheng et al. [26] investigated influence of the thermal boundary conditions on the three-dimensional convective flow in a vertical cylinder heated from below. These authors revealed that in case of the adiabatic side surface the flow is highly asymmetric and contains multicellular vortices even at symmetric mathematical statement of the problem.
Fig. 2. Variation of the average Nusselt number versus the dimensionless time and the mesh parameters at Ra = 104, Pr = 0.7, k2,1 = 5.7 104, c = p/6.
3584
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 3. Velocity field and isosurface of the vertical velocity component ~ V and temperature field H at k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 for Ra = 104 (a), Ra = 5 104 (b) and Ra = 105 (c).
The objective of the current investigation is mathematical simulation of transient natural convection modes in an inclined cylinder having heat-conducting walls of finite thickness at presence of a local heat source in conditions of convective heat exchange with an environment.
2. Physical and mathematical model The physical model of transient conjugate natural convection in an inclined closed cylinder and the coordinate system are schemat-
ically shown in Fig. 1. The domain of interest consists of a fluid cavity bounded by heat-conducting walls of finite thickness. A heat source located at the bottom of the cavity is kept at constant temperature. The convective heat exchange with an environment is modeled at borders r = Lr and z = z1. The border z = 0 is supposed to be adiabatic. It is assumed in the analysis that the thermophysical properties of the walls material and of the fluid are independent of temperature, and the flow is laminar. The fluid is viscous, heat-conducting, Newtonian, and the Boussinesq approximation is valid. The fluid motion and heat transfer in the cavity are supposed to be
3585
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
three-dimensional. Radiation heat exchange from the heat source and between the walls is neglected in comparison with convection. Taking into account the above-mentioned assumptions, the differential problem of non-stationary conjugate natural convection in an inclined closed cylinder on the basis of conservation laws of mass, momentum and energy can be written as follows [27]: In the fluid cavity
aim of the present investigation is the analysis of a thermal state of the system in conditions of a mutual influence of natural convection inside a vertical cylinder filled with fluid and heat conduction in solid walls. For this purpose we shall enter into consideration satisfying the condition ¼ w i þ w j þ w k, the vector potential w r u z ¼ 0 [29,30]. Then components of the velocity will such as div w be defined as follows:
@ðrV r Þ @V u @ðrV z Þ þ þ ¼ 0; @r @z @u
¼ rot w V 1 @wz @wu @w @w 1 @ðrwu Þ 1 @wr ) Vr ¼ : ; Vu ¼ r z ; Vz ¼ r @u r @r r @u @z @z @r
ð1Þ
2 @V r @V r V u @V r @V r V u þ Vz þ Vr þ @t @r r @u @z r @V 1 @p V 2 u r gbðT T 0 Þ þ m r2 V r 2 2 ¼ r @u q2 @r r
sin c cos u;
ð7Þ The modified vorticity functions will become
ð2Þ
xr ¼ r
@V u @V z ; @z @u
xu ¼
@V z @V r ; @r @z
xz ¼
@V u @V r Vu r : @u @r ð8Þ
@V u @V u V u @V u @V u V r V u þ Vz þ Vr þ þ @t @r r @u @z r V u 2 @V r 1 @p 2 þ gbðT T 0 Þ sin c sin u; þ m r Vu 2 þ 2 ¼ r q2 @ u r @u r ð3Þ @V z @V z V u @V z @V z þ Vz þ Vr þ @t @r r @u @z 1 @p þ mr2 V z þ gbðT T 0 Þ cos c; ¼ q2 @z @T @T V u @T @T þ Vr þ ¼ a2 r2 T; þ Vz @t @r @z r @u
ð4Þ
Z ¼ z=Lr ; s ¼ tV b =Lr ; H ¼ ðT T 0 Þ=ðT hs T 0 Þ; 8 8 8 > < U ¼ V r =V b ; > < Wr ¼ wr =ðV b Lr Þ; > < Xr ¼ xr =V b ; 2 2 2 ~ V ¼ V u =V b ; Wu ¼ wu =ðV b Lr Þ; Xu ¼ xu Lr =V b ; r ¼ Lr r ; > > > : : : W ¼ V z =V b ; Wz ¼ wz =ðV b Lr Þ; Xz ¼ xz =V b :
ð5Þ
The dimensionless equations for the fluid can be written as follows:
R ¼ r=Lr ;
~ 2 Wr Wr 2 R2 r
In the solid walls
@T ¼ a1 r2 T: @t
These functions can be derived by means of the vorticity vector ¼ rot V ~ ¼x ~ uj þ x ~ zk ~ ri þ x like that components x ~ r ; xu ¼ x ~ u ; xz ¼ r x ~ z. xr ¼ rx A mathematical model was formulated in terms of the dimensionless variables the vector potential and the modified vorticity functions by using the non-dimensional variables:
ð6Þ
A transformation of the formulated system of differential equations (1)–(6) to a form eliminating direct search of the pressure field [21,28] is represented as the most reasonable because the
@ Wu ¼ RXr ; @u
~ 2 Wu Wu þ 2 R2 r ~ 2 Wz ¼ Xz ; Rr
@ Wr ¼ R2 Xu ; @u
ð9Þ
ð10Þ ð11Þ
Fig. 4. Isolines of vector potential component Wr and temperature fields H at R = 0.28, k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 for Ra = 104 (a), Ra = 5 104 (b) and Ra = 105 (c).
3586
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 5. Isolines of vector potential component Wu and temperature fields H at u = 0.56p, k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 for Ra = 104 (a), Ra = 5 104 (b) and Ra = 105 (c).
@ Xr @ Xr V @ Xr @ Xr @U @U @U U Xr @V Xu þU þW Xz þ Xr U @R @u @Z @Z @s @R R @ u @Z R rffiffiffiffiffiffi @ X Pr ~ 2 2 @ Xr 2 @H @H u þR r Xr cos c; ¼ sin c sin u Ra R @R R @ u @Z @u ð12Þ
@ Xu @ Xu V @ Xu @ Xu Xr @V Xu @V Xz @V U Xu V Xr þ þ 2 þU þW R @u @s @R @Z R @R R @ u R @Z R R rffiffiffiffiffiffi Xu 2 @ Xr Pr ~ 2 @H @H ¼ r Xu 2 þ 3 þ sin c cos u þ cos c; @Z @R Ra R @u R ð13Þ
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3587
Fig. 6. Isolines of vector potential component Wz and temperature fields H at Z = 0.56, k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 for Ra = 104 (a), Ra = 5 104 (b), Ra = 105 (c).
@ Xz @ Xz V @ Xz @ Xz @W @W @W U Xz @V þU þW Xz þ Xr þU Xu @s @R R @ u @Z R @R @u @Z @R rffiffiffiffiffiffi Pr ~ 2 X 2 @ Xz @H @H r Xz þ 2z sin c cos u R ¼ sin c sin u; @u @R Ra R @R R ð14Þ
@H @H V @H @H 1 ~ 2 H: þU þW þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi r @s @R R @ u @Z Pr Ra
ð15Þ
Energy equation for solid walls:
1 @H ~ 2 ¼ r H: Fo @ s
ð16Þ
Eqs. (9)–(11) have been obtained as a result of substitution of the velocity components (7) in the definition of the modified vorticity functions (8). Eqs. (12)–(14) have been obtained according to the entered modified vorticity functions (8). For example, Eq. (12) has been obtained as follows. Firstly we differentiate Eq. (3) with
3588
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Initial conditions are
Wr ðR; u; Z; 0Þ ¼ 0; Wu ðR; u; Z; 0Þ ¼ 0; Wz ðR; u; Z; 0Þ ¼ 0; Xr ðR; u; Z; 0Þ ¼ 0; Xu ðR; u; Z; 0Þ ¼ 0; Xz ðR; u; Z; 0Þ ¼ 0; H(R, u, Z, 0) = 0 except the temperature for the heat source at which H = 1 during the whole process time. Boundary conditions are: convective heat exchange with an environment is modeled at the walls such as R = 1 and Z = z1/Lr: @@Hn ¼ BiðHe HÞ; at the wall Z = 0 for Eq. (16) the heat insulation condition is set
@H ¼ 0; @Z at the solid–fluid interface R = r1/Lr:
@ Wr ¼ Wu ¼ Wz ¼ 0; @R
(
H1 ¼ H2 ; @ H1 @R
H2 ¼ k2;1 @@R ;
at the solid–fluid interfaces Z = const except the heat source surface: Fig. 7. Variation of the average Nusselt number versus the Rayleigh number and the inclination angle at k2,1 = 5.7 104, s = 100, Pr = 0.7
respect to z. The received equation is multiplied by r and Eq. (4) differentiated with respect to u is subtracted from it. Eqs. (9)–(16) are subjected to the following initial and boundary conditions.
@ Wz Wr ¼ Wu ¼ ¼ 0; @Z
(
H1 ¼ H2 ; @ H1 @Z
H2 ¼ k2;1 @@Z ;
Wz at the heat source surface: Wr ¼ Wu ¼ @@Z ¼ 0; H ¼ 1:0; hs for the heat source ðat lLwr 6 Z 6 lw þz ; 0 6 R 6 rL1r ; 0 6 u < 2pÞ Lr H = 1.0; at u = 0, 2p the periodicity conditions are set, e.g.
Fig. 8. Isolines of vector potential component Wr and temperature fields H at R = 0.28, k2,1 = 5.7 104, Ra = 104, Pr = 0.7, c = p/3 for s = 3 (a), s = 6 (b), s = 20 (c) and s = 70 (d).
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
8 8 8 < Wr ju¼0 ¼ Wr ju¼2p ; < Xr ju¼0 ¼ Xr ju¼2p ; < Hju¼0 ¼ Hju¼2p ; @ H : @@Wur ¼ @@Wur ; : @@Xur ¼ @@Xur ; : @@H ¼ : u @ u u¼0
u¼2p
u¼0
u¼2p
u¼0
u¼2p
The correct boundary conditions at R = 0 for vector potential components and modified vorticity functions and also for temperature have been formulated in accordance with the differential equations. The detailed description of such boundary conditions has been presented by Paasonen [31]. Boundary conditions for the modified vorticity functions are obtained on the basis of the following correlations:
Xr ¼ R
@V @W ; @Z @ u
Xu ¼
@W @U ; @R @Z
Xz ¼
@U @V : V R @u @R
Converting these correlations for the internal boundaries such as R = const, Z = const one can obtain: – Xr ¼ 0; Xu ¼ @W ; Xz ¼ R @V at R ¼ const; @R @R @U – Xr ¼ R @V ; X ¼ ; X ¼ 0 at Z ¼ const: u z @Z @Z
3589
Thus for definition of the modified vorticity functions at internal boundaries it is necessary to determine the velocity field: @W U ¼ 1R @@Wuz @Zu ;
V¼
@ Wr @ Wz ; @Z @R
W¼
1 @ðRWu Þ 1 @ Wr : R @R R @u
Eqs. (9)–(16) with corresponding initial and boundary conditions has been solved by means of finite differences method [21,28,32,33]. Uniform grid in R, u and Z directions was used for all computations. A locally one-dimensional scheme of Samarskii [32] was used to solve Eqs. (12)–(16) numerically. In this scheme a solution of the three-dimensional scheme reduces to a sequential solution of the one-dimensional systems. It is to be noted here that Samarskii monotone approximation was applied to discretization of the convective terms in Eqs. (12)–(15) [33]. An implicit difference scheme was used. An evolutionary term represented a one-sided difference in time and had the first order of accuracy in time step. All derivatives with respect to spatial coordinates were approximated with
Fig. 9. Isolines of vector potential component Wu and temperature fields H at u ¼ 0:56p, k2;1 ¼ 5:7 104 ; Ra ¼ 104 ; Pr ¼ 0:7; c ¼ p=3 for s = 3 (a), s = 6 (b), s = 20 (c) and s = 70 (d).
3590
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 10. Isolines of vector potential component Wz and temperature fields H at Z ¼ 0:56, k2;1 ¼ 5:7 104 ; Ra ¼ 104 ; Pr ¼ 0:7; c ¼ p=3 for s = 3 (a), s = 6 (b), s = 20 (c) and s = 70 (d).
the second order of accuracy in the step along the coordinate. As a result an error of approximation of the partial differential equations and the boundary conditions was O(Ds + (Dr)2 + (D
u)2 + (Dz)2). A computational time step was selected from conditions of a stable algorithm and so Ds = 103 [21,28,34]. It should be noted that analysis of the time step effect was previously
3591
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
8 W W W Wzi;j1;k > U i;j;k ¼ zi;jþ1;k ui;j;kþ12hz ui;j;k1 ; > 2Ri hu > < W W W W V i;j;k ¼ ri;j;kþ12hz ri;j;k1 ziþ1;j;k2hr zi1;j;k ; > > > : W ¼ Wuiþ1;j;k Wui1;j;k þ Wui;j;k Wri;jþ1;k Wri;j1;k : i;j;k 2hr R 2R hu i
Fig. 11. Variation of the average Nusselt number versus the dimensionless time and the Rayleigh number at Pr = 0.7, k2,1 = 5.7 104, c = p/3.
performed with time steps Ds of 0.01, 0.001, 0.0001. It was found that at Ds = 102 increasing oscillations of the average Nusselt number have been observed at high analyzed Rayleigh numbers. The numerical stability has been marked for time steps of 0.001 and 0.0001. Therefore the time step of 0.001 was adopted in the present study. Eqs. (9)–(11) have been discretised by means of seven-points difference scheme, so-called ‘‘cross scheme’’, on the basis of a symmetric approximation of the second derivatives. The obtained system of linear algebraic equations has been solved by the successive over relaxation method. Optimum value of a relaxation parameter has been chosen on the basis of computing experiments. Vorticity boundary conditions have an essential value for numerical simulation of analyzed process. Approximation of these boundary conditions for the internal boundaries such as R = r1/Lr, Z = (lw + zhs)/Lr, Z = (z1 - lw)/Lr with the second order of accuracy on the basis of expansion in a Taylor series can be written as follows
8 Xrnr;j;k ¼ 0; > > < 3W 4W nr1;j;k þW nr2;j;k Xunr;j;k ¼ @W ¼ nr;j;k ; at R ¼ Rnr ¼ r 1 =Lr ; — @R nr;j;k 2hr > > 4V nr1;j;k V nr2;j;k 3V nr;j;k : @V Xznr;j;k ¼ R @R nr;j;k ¼ Rnr 2hr 8 4V i;j;nz þ1 V i;j;nz þ2 3V i;j;nz > 1 1 1 > X ¼ R @V ¼ Ri ; > 2hz @Z i;j;nz1 < ri;j;nz1 lw þ zhs 3U þU 4U i;j;nz1 i;j;nz1 þ2 i;j;nz1 þ1 — @U X ¼ @Z i;j;nz ¼ ; at Z ¼ Z nz1 ¼ Lr ; > 2hz 1 > ui;j;nz1 > : Xzi;j;nz1 ¼ 0 8 @V 3V i;j;nz 4V i;j;nz 1 þV i;j;nz 2 > 2 2 2 > ; > Xri;j;nz2 ¼ R @Z i;j;nz2 ¼ Ri 2hz < z 1 lw 4U U 3U i;j;nz 1 i;j;nz 2 i;j;nz — 2 2 2 Xui;j;nz2 ¼ @U ¼ ; atZ ¼ Z nz2 ¼ Lr ; > 2h @Z i;j;nz z 2 > > : Xzi;j;nz2 ¼ 0 where i, j, k are mesh point coordinates. For definition of the velocity field the second-order central differences relations are used for spatial derivatives such as
i
The concise description of numerical simulation of boundary conditions at a cylinder axis in the present work is impossible because the numerical approximation becomes complicated not only owing to the analysis of a nonsymmetric case, but also owing to the selected variables such as vector potential components and modified vorticity functions. The detailed description of the chosen method has been completely presented by Paasonen [31]. It should be noted that for realization of the coordinate-splitting method (a locally one-dimensional scheme of Samarskii [32]) the correct sequence in a choice of coordinates is important. Firstly, the solution along circles with the periodicity conditions at u = 0 and u = 2p is carried out, then the solution along all internal grid lines parallel to an axis of the cylinder, except for this axis, with conditions at the bases of cylinder is conducted. Further on at the cylinder axis the equation obtained by means of the coordinate-splitting method is solved. At the third stage there is a solution of the system of the one-dimensional equations along the radiuses of the cylinder by means of the method presented in [31]. The method of solution has been tested on different meshes. Fig. 2 shows sensitivity of the average Nusselt number at the heat R 0:88 R 2p @ H 1 source surface Nuavg ¼ 1:76 j @Z jZ¼0:24 du dR to the mesh 0 p 0 parameters with increase in time at Ra = 104, Pr = 0.7, c = p/6, k2,1 = 5.7 104. On the ground of optimization of calculation accuracy for 0 6 s 6 100 the computational mesh such as 50 50 50 has been selected. There were 10 nodes along the R-coordinate inside the solid walls. 3. Results and discussion Numerical investigation of the boundary value problem (9)– (16) has been carried out at following values of dimensionless complexes: Ra = 104, 5 104, 105; Pr = 0.7, 7.0; 0 6 s 6 100; c = 0, p/6, p/3, p/2; He = 1; k2,1 = 5.7 104 – Bi = 4.35 102, Fo = 3.5 104; and k2,1 = 4.3 102 – Bi = 3.33, Fo = 9.38 106. Analysis has been carried out for the following geometrical parameters: z1/Lr = 2, lw/Lr = 0.12, r1/Lr = 0.88, zhs/Lr = 0.12. Particular efforts have been focused on the effects of four types of influential factors such as the Rayleigh and Prandtl numbers, the thermal conductivity ratio and an inclination angle on the thermohydrodynamic regimes. 3.1. Effect of Rayleigh number Velocity and temperature fields and isosurfaces of the vertical velocity component at k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 and different values of the Rayleigh number are presented in Fig. 3. This figure reflects formation of the vortex structures in the gas cavity having the inclination angle c = p/3. The motion paths of gas particles crossing the cylinder axis demonstrate the presence of inclination angle. More detailed analysis of distributions of thermohydrodynamic parameters will be presented for specific sections of the domain of interest. Figs. 4–6 show the isolines of the vector potential components and the isotherms for various cross-sections at k2,1 = 5.7 104, s = 100, Pr = 0.7, c = p/3 and different values of the Rayleigh number Ra = 104, 5 104, 105. The increase in the Rayleigh number for the considered range is reflected in an insignificant modification of the vector potential component Wr isolines (Fig. 4). Two convective cells are formed
3592
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 12. Velocity field and isosurfaces of the vertical velocity component ~ V and temperature field H at Ra = 104, k2,1 = 5.7 104, s = 100, c = p/6 for Pr = 0.7 (a) and Pr = 7.0 (b).
in the cross-section R = 0.28. These cells differ in the direction of fluid motion and reflect formation of a thermal plume symmetrically with respect to axis u = p. It is necessary to note the displacement of a convective cells core to the axis u = p with the increase in the Rayleigh number that is caused by an intensification of the transport processes of momentum and energy in this zone. Appearance of a conductive heat transfer zone near the heat source surface is obvious. It is confirmed by the isotherms position. The increase in the Rayleigh number also leads to a decrease in the qffiffiffiffi
Pr thermal boundary layer thickness O 4 Ra [35] above the heat source. In all probability, the formation of a thermal plume and convective cells symmetrically with respect to the axis u = p, is caused by the presence of the cylinder inclination angle such as c = p/3. Interaction of the low temperature front spreading from the border Z = z1/Lr deep into the domain of interest, with the
thermal plume is observed also. Thus an increase in the temperature differences, owing to a rise in a heat source temperature, leads to a decrease in the sizes of the zone that is subjected to the influence of the low temperature. The increase in the Rayleigh number results not only in an increment of the temperature for the upper layers of the analyzed object, but also in an increase in a thermal plume thickness that is reflected in formation of a significant cross-section energy flow. Analyzing changes of the vector potential component Wu for the cross-section u = 0.56p (Fig. 5) it is possible to note an increase in sizes and intensity of a circulating flow located in the angular zone of the domain of interest. The latter leads to distortion of the basic convective cell and to reduction of the fluid circulation velocities in it. Displacement of a core of the basic vortex is caused by both a development of a recirculation flow described above and a formation of a thermal plume over the heat source close to the solid walls.
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3593
Fig. 13. Isolines of vector potential component Wr and temperature fields H at R = 0.28, k2,1 = 5.7 104, Ra = 104, c = p/6, s = 100 for Pr = 0.7 (a) and Pr = 7.0 (b).
Distribution of the isotherms in cross-section u = 0.56p as well as at R = 0.28 reflects both an increase in the temperature for the upper layers of the domain of interest that interferes with penetration of the low temperature front and a drift of convective cell core to zone 0.5 < R < 0.7. The presence of thermal spots in the cavity characterizes an essential three-dimensionality of the analyzed process. Distributions of the vector potential component Wz and the isotherms for the cross-section Z = 0.56 (Fig. 6) confirm the presence of symmetry with respect to the axis u = p observable in Fig. 4, and also reflect both an intensification of convective flow at the increase in the Rayleigh number and a heating of the fluid cavity. The analysis of the influence of the Rayleigh number and the inclination angle on the average Nusselt number at the heat source surface (Fig. 7) has been conducted. An increase in the generalized heat transfer coefficient at the increment in the Rayleigh number is observed. The increase in the inclination angle leads to convective heat transfer intensification owing to an increase in the temperature gradient at the heat source surface. The nonmonotonic changes in the average Nusselt number at a variation of the inclination angle a reflects the presence of an optimum position of the cylinder with the maximum value of Nuavg. In the analyzed range of the Rayleigh number it is possible to indicate two intervals with the maximum magnitude of the generalized heat transfer coefficient at various values of the inclination angle of the cylinder such as at 104 6 Ra 6 3 104 we have c = p/3 and at 3 104 < Ra 6 105 we have c = p/6. The presented average Nusselt number curves in Fig. 7 have been defined as correlations that are very convenient in carrying out applied calculations:
Nuavg ¼ 0:188 Ra0:315 at c ¼ 0; Nuavg ¼ 0:242 Ra0:305 at c ¼ p=6; Nuavg ¼ 0:323 Ra0:277 at c ¼ p=3; Nuavg ¼ 0:34 Ra0:265 at c ¼ p=2:
3.2. Effect of the transient factor
Fig. 14. Isolines of vector potential component Wu and temperature fields H at u = 0.56p, k2,1 = 5.7 104, Ra = 104, c = p/6, s = 100 for Pr = 0.7 (a) and Pr = 7.0 (b).
The transient factor in conjugate heat transfer problems plays an essential role [21,28] as it defines not only dynamics of the velocity and temperature fields owing to formation, evolution and dissipation of the vortex structures in a fluid cavity, but also a thermal lag of solid walls in conditions of heat exchange with an environment.
3594
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 15. Isolines of vector potential component Wz and temperature fields H at Z = 0.56, k2,1 = 5.7 104, Ra = 104, c = p/6, s = 100 for Pr = 0.7 (a) and Pr = 7.0 (b).
Fig. 8 demonstrates isolines of the vector potential component
Wr and temperature fields H in the cross-section R = 0.28 at k2,1 = 5.7 104, Ra = 104, Pr = 0.7, c = p/3 for different values of dimensionless time. Two convective cells of small intensity are formed symmetrically with respect to axis u = p at the initial time s = 3 in cross-section R = 0.28 that is caused by conductive heat transfer in the cavity (isotherms are parallel to the heat source surface). The increase in dimensionless time leads to formation and evolution of a thermal plume which intensifies convective heat transfer in the cavity. The sizes of convective cells are increased while cores of the vortexes are displaced to a center layer of the cylinder. Penetration of a low temperature front from the border Z = z1/Lr deep into the cavity is observed jointly with evolution of the thermal plume. In the cross-section u = 0.56p (Fig. 9) the dynamics of formation and evolution of the hydrodynamic structures also reflects the basic moments described above. The increase in the dimensionless time s leads to an intensification of the basic convective cell formed directly above the heat source. An appearance of a vortex is observed in the zone 1 < Z < 1.8 owing to cooling of the top wall of the domain of interest. This vortex loses an internal energy and decreases in sizes at the increase in dimensionless time. Distribution of isotherms in the fluid cavity at s P 20 is defined not only by influence of the buoyancy force, but also by dynamics of con-
Fig. 16. Variation of the average Nusselt number versus the dimensionless time and the Prandtl number at Ra = 104, k2,1 = 5.7 104, c = p/6.
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
ductive heat transfer in the solid walls of the cylinder. Presence of the inclination angle such as c = p/3 and heating of the solid walls lead to formation of the thermal plume near to the side surface of the cylinder. The modification of the isolines of the vector potential component Wz in the cross-section Z = 0.56 (Fig. 10) reflects an essential nonstationarity of the analyzed process. Formation of two groups of low-intensity convective cells differing by sizes is observed at the initial time. Disintegration of the large vortexes to the finer vortexes and reorientation of formed hydrodynamic structures are observed at the increase in the dimensionless time. There is an increase in sizes of the central vortexes group at s = 20; and by the moment s = 70 this two-vortical system becomes dominating. It is necessary to note, that the hydrodynamic
3595
structures formed at the moments s = 3 and s = 70 are similar to each other. A basic difference consists in orientation and intensity of these vortexes. The temperature field varies insignificantly. An increase in the thickness of thermal plume located near to the side surface of the cylinder in the zone u = p is observed. Dynamics of the generalized heat transfer coefficient at the heat source surface at various values of the Rayleigh number is presented in Fig. 11. Change of the average Nusselt number can be broken into two levels. 0 < s < 20 is an initial level reflecting presence of oscillations in a function of Nuavg(s). 20 < s < 100 is a level of a thermal stabilization of natural convection. It is necessary to note, that at large Rayleigh numbers ðRa P 105 Þ stabilization is delayed owing to essential interaction of the temperature fields of solid walls and the fluid cavity.
Fig. 17. Velocity field and isosurfaces of the vertical velocity component ~ V and temperature field H at Ra = 105, Pr = 0.7, c = p/6, s = 100 for k2,1 = 5.7 104 (a) and k2,1 = 4.3 102 (b).
3596
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3.3. Effect of the Prandtl number Further the effect of fluid heat-transfer properties on the flow and heat transfer modes will be analyzed. The increase in the Prandtl number in 10 times leads to attenuation of convective heat transfer in the cavity (Figs. 12–15,b), that is caused by an increase in a kinematic viscosity coefficient and an internal friction force. The latter is reflected as well in both velocity field (Fig. 12) and configurations of isolines of the vector potential components e.g. Wz (Fig. 15). Slight displacement of the convective cells cores at Pr = 7.0 in the cross-section R = 0.28 (Fig. 13) is caused by more intensive diffusion of perturbations in the cavity from solid walls. An increase in sizes of a recirculation zone located in an angular area of the fluid cavity is observed in the cross-section u = 0.56p (Fig. 14). The latter is related also to heat conduction dynamics in solid walls. It is possible to note more intensive heating of the domain of interest at Pr = 7.0 (Fig. 12) that is caused by both hydrodynamics, namely, small flow velocities, presence of recirculation zones reflecting formation of high temperature areas, and thermodynamics of the conjugate problems, namely, interaction of conductive heat transfer in solid walls and natural convection in the fluid cavity. The latter leads to a decrease in an intensity of the low temperature penetration from the border Z = z1/Lr deep into the domain of interest (Figs. 12 and 14), and also to an increase in a thermal plume thickness (Figs. 12–15). It is necessary to note formation of the thermal spots (Figs. 13 and 14), that is caused by the three-dimensionality of the analyzed process. Change dynamics of the generalized heat transfer coefficient versus the Prandtl number is presented in Fig. 16. An increase in the fluid viscosity is reflected in a delaying of the thermal stabilization zone. At initial time s < 47 a maximum intensity of convective heat transfer is observed at Pr = 7.0, and at s > 47 vice versa. 3.4. Effect of the thermal conductivity ratio Presence of heat-conducting walls of finite thickness can essentially modify the temperature field in the fluid cavity and flow regimes [12,15]. The solid walls can play a role as an original buffer zone owing to a thermal lag of the solid material. For a concrete time interval walls can accumulate energy, reflecting indirect influence of an environment on the processes proceeding in the cavity. Figs. 17–20 show the velocity and temperature fields and also the isolines of the vector potential components in various crosssections, corresponding to convective heat transfer mode
Fig. 19. Isolines of vector potential component Wu and temperature fields H at u = 0.56p, Ra = 105, Pr = 0.7, c = p/6, s = 100 for k2,1 = 5.7 104 (a) and k2,1 = 4.3 102 (b).
Fig. 18. Isolines of vector potential component Wr and temperature fields H at R = 0.28, Ra = 105, Pr = 0.7, c = p/6, s = 100 for k2,1 = 5.7 104 (a) and k2,1 = 4.3 102 (b).
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3597
Fig. 20. Isolines of vector potential component Wz and temperature fields H at Z = 0.56, Ra = 105, Pr = 0.7, c = p/6, s = 100 for k2,1 = 5.7 104 (a) and k2,1 = 4.3 102 (b).
Ra = 105, Pr = 0.7, s = 100, c = p/6, at different values of the thermal conductivity ratio k2,1 = 5.7 104, 4.3 102. The increase in the thermal conductivity ratio is caused with a decrease in the thermal conductivity of the solid material. Transition from k2,1 = 5.7 104 to k2,1 = 4.3 102 leads to essential modification of the vector potential component Wz in the cross-section Z = 0.56 (Fig. 20). Combination of the vortical structures with formation of the global circulation defining mass transfer, momentum transport and energy transport in the cavity is observed. The flow intensity is increased, though slight decrease in the fluid velocities is observed in a vertical direction (Fig. 18). The most essential changes concern the temperature fields. The temperature in the fluid cavity is decreased at k2,1 = 4.3 102 (Fig. 17) that is related to a diminution of the thermal conductivity of solid walls and, accordingly, with low speed of the energy transport in solid walls both from the heat source and from an environment. The latter is evidently presented in Fig. 19 where the isotherm corresponding to the dimensionless temperature H = 0.2 at k2,1 = 4.3 102 is inside of the top wall and in case of k2,1 = 5.7 104 it is at the solid–fluid interface. It is necessary to note loss of stability of the thermal plume above the energy source (Fig. 19). It reflects formation of transition thermohydrodynamic regimes. More intensive cooling of the upper wall in case of a heat-conductive material leads to an increase in a temperature gradient at the internal solid–fluid interface. The latter leads to an increase in convective flow intensity in the cavity.
Fig. 21. Variation of the average Nusselt number versus the Rayleigh number and the thermal conductivity ratio at Pr = 0.7, c = p/6, s = 100.
3598
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
Fig. 22. Isolines of vector potential component Wr and temperature fields H at R = 0.28, Ra = 104, Pr = 0.7, k2,1 = 5.7 104, s = 100 for c = 0 (a) and c = p/2 (b).
Symmetry of the vector potential components isolines and isotherms with respect to an axis u = p is conserved (Fig. 18), that is caused by both the inclination angle of the cylinder and with indirect influence of an environment. Variation of the average Nusselt number at the heat source surface depending on the Rayleigh number in a range 104 6 Ra 6 105 and the thermal conductivity ratio is presented in Fig. 21. The increase in k2,1 leads to formation of two zones reflecting various behavior of the generalized heat transfer coefficient. The increase in the thermal conductivity ratio leads to an increase in the average Nusselt number at 104 6 Ra < 5:6 104 owing to rising of the temperature gradient that is caused by the growth of energy removal in solid walls. The opposite situation caused by an increase in the role of the hydrodynamic component of the analyzed process is observed in case of 5:6 104 < Ra 6 105 . An increase in the buoyancy force leads to an intensive vertical energy transport in the cavity i.e. the near-surface layers close to the heat source are considerably heated. 3.5. Effect of the inclination angle
Fig. 23. Isolines of vector potential component Wu and temperature fields H at u = 0.56p, Ra = 104, Pr = 0.7, k2,1 = 5.7 104, s = 100 for c = 0 (a) and c = p/2 (b).
Figs. 22–24 present the isolines of vector potential components and the temperature fields in the various cross-sections, corresponding to the convective heat transfer regime at Ra = 104, Pr = 0.7, k2,1 = 5.7 104, s = 100 for various values of the inclination angle c = 0 [34] and c = p/2 . Presence of the inclination angle such as c = p/2 leads to displacement of a thermal plume to an axis u = p and in an intensification of convective flow in the cavity (Fig. 22). Circulating flows reflecting a fluid motion in an opposite direction in comparison with the regime at c = 0 [34] are shaped in the cross-section u = 0.56p (Fig. 23). Clear differences in hydrodynamic structures are displayed in the cross-section Z = 0.56 (Fig. 24) reflecting a system rotation by the angle c = p/2. The temperature fields also differ. There is more intensive cooling of the cavity from an environment in the crosssection u = 0.56p at c = 0. Variations of the average Nusselt number versus the dimensionless time at different values of the inclination angle are presented in Fig. 25. It is necessary to note on the one hand an essential increase in Nuavg at presence of an inclination angle of the cylinder and on the other hand a smoother stabilization of the process. While reorganization of the thermodynamic structures at the moment s = 58 at c = 0 defining presence of a breakpoint in dependence Nuavg(s) is observed. The increase in an inclination
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
3599
Fig. 24. Isolines of vector potential component Wz and temperature fields H at Z = 0.56, Ra = 104, Pr = 0.7, k2,1 = 5.7 104, s = 100 for c = 0 (a), c = p/2 (b).
angle leads to an increase in the generalized heat transfer coefficient in a zone s < 10, subsequently this dynamics is broken, and also is reflected in faster stabilization of the analyzed process. 4. Conclusions
Fig. 25. Variation of the average Nusselt number versus the dimensionless time and the inclination angle of the cylinder at Ra = 104, Pr = 0.7, k2,1 = 5.7 104.
Transient three-dimensional natural convection in an inclined closed tube with heat-conducting walls of finite thickness at presence of the internal heat source of the constant temperature in conditions of convective heat exchange with an environment has been numerically solved. Distributions of isolines of the vector potential components and the temperature fields describing features of the analyzed flow and heat transfer regimes at Ra = 104, 5 104, 105, Pr = 0.7, 7.0, 0 6 s 6 100, c = 0, p/6, p/3, p/2, k2,1 = 5.7 104, 4.3 102 have been obtained. A detailed analysis of influence of the Rayleigh number and the Prandtl number, of the dimensionless time, of the inclination angle of the tube and of the thermal conductivity ratio both on thermohydrodynamic fields and on the average Nusselt number at the heat source surface has been conducted. It has been ascertained that in the analyzed range of the Rayleigh number it is possible to indicate two intervals with a maximum magnitude of the generalized heat transfer coefficient at various values of the inclination angle of the tube such as at 104 6 Ra 6 3 104 we have c = p/3 and at 3 104 < Ra 6 105 we have c = p/6 (Fig. 7). It has been shown that variation of the average Nusselt number can be broken into two levels (Fig. 11).
3600
M.A. Sheremet / International Journal of Heat and Mass Transfer 55 (2012) 3582–3600
0 < s < 20 is an initial level reflecting presence of oscillations in a function Nuavg(s). 20 < s < 100 is a level of thermal stabilization of natural convection. The increase in the Prandtl number is reflected in a delaying of the thermal stabilization zone (Fig. 16). At initial time s < 47 maximum intensity of convective heat transfer is observed at Pr = 7.0, and at s > 47 maximum intensity of convective heat transfer is observed at Pr = 0.7. The various mechanisms defining change of the average Nusselt number versus the thermal conductivity ratio have been ascertained (Fig. 21). The increase in k2,1 leads to an increase in Nuavg at 104 6 Ra < 5:6 104 owing to a prevailing influence of thermal component of the analyzed process. An opposite situation caused by an increase in a role of the hydrodynamic component of natural convection is observed in case of 5:6 104 < Ra 6 105 . An essential increase in Nuavg and a smoother stabilization of the process at presence of the inclination angle of the cylinder have been shown (Fig. 25). Acknowledgements This work was performed according to the Russian Federal Purpose Program during 2009–2013 years (state contract No. P 357) and was supported by the Grants Council (under the President of the Russian Federation), Grant No. MR-5652.2012.8. References [1] G. Carbajal, C.B. Sobhan, G.P. ‘‘Bud’’ Peterson, D.T. Queheillalt, Haydn N.G. Wadley, A quasi-3D analysis of the thermal performance of a flat heat pipe, Int. J. Heat Mass Transfer 50 (2007) 4286–4296. [2] L.L. Vasiliev, Heat pipes in modern heat exchangers, Appl. Therm. Eng. 25 (2005) 1–19. [3] G.V. Kuznetsov, M.A. Al-Ani, M.A. Sheremet, Numerical analysis of convective heat transfer in a closed two-phase thermosyphon, J. Eng. Thermophys. 20 (2011) 201–210. [4] Gilles Desrayaud, Alberto Fichera, Manuel Marcoux, Numerical investigation of natural circulation in a 2D-annular closed-loop thermosyphon, Int. J. Heat Fluid Flow 27 (2006) 154–166. [5] L.A. Florio, A. Harnoy, Combination technique for improving natural convection cooling in electronics, Int. J. Therm. Sci. 46 (2007) 76–92. [6] H. Oprins, J. Danneels, B. Van Ham, B. Vandevelde, M. Baelmans, Convection heat transfer in electrostatic actuated liquid droplets for electronics cooling, Microelectron. J. 39 (2008) 966–974. [7] A.N. Campbell, S.S.S. Cardoso, A.N. Hayhurst, A comparison of measured temperatures with those calculated numerically and analytically for an exothermic chemical reaction inside a spherical batch reactor with natural convection, Chem. Eng. Sci. 62 (2007) 3068–3082. [8] Y. Jaluria, Design and Optimization of Thermal Systems, McGraw-Hill, New York, 1998. p. 626. [9] V.I. Polezhaev, Free convection in the internal problem: results and prospects, J. Eng. Phys. Thermophys. 69 (1996) 676–687. [10] M.V. Diomidov, M.I. Nizovtsev, V.I. Terekhov, Ventilation of window interplane cavity aimed at a higher temperature of the inner pane, Therm. Sci. 6 (2002) 15–22. [11] D.A. Kaminski, C. Prakash, Conjugate natural convection in a square enclosure effect of conduction on one of the vertical walls, Int. J. Heat Mass Transfer 29 (1986) 1979–1988.
[12] A. Liaqat, A.C. Baytas, Numerical comparison of conjugate and non-conjugate natural convection for internally heated semi-circular pools, Int. J. Heat Fluid Flow 22 (2001) 650–656. [13] A. Liaqat, A.C. Baytas, Conjugate natural convection in a square enclosure containing volumetric sources, Int. J. Heat Mass Transfer 44 (2001) 3273– 3280. [14] G.V. Kuznetsov, M.A. Sheremet, Conjugate heat transfer in an enclosure under the condition of internal mass transfer and in the presence of the local heat source, Int. J. Heat Mass Transfer 52 (2009) 1–8. [15] M.A. Sheremet, The influence of cross effects on the characteristics of heat and mass transfer in the conditions of conjugate natural convection, J. Eng. Thermophys. 19 (2010) 119–127. [16] Yasin Varol, Hakan F. Oztop, Ahmet Koca, Filiz Ozgen, Natural convection and fluid flow in inclined enclosure with a corner heater, Appl. Therm. Eng. 29 (2009) 340–350. [17] Chengwang Lei, S.W. Armfield, J.C. Patterson, Unsteady natural convection in a water-filled isosceles triangular enclosure heated from below, Int. J. Heat Mass Transfer 51 (2008) 2637–2650. [18] Fu-Yun Zhao, Di Liu, Guang-Fa Tang, Natural convection in an enclosure with localized heating and salting from below, Int. J. Heat Mass Transfer 51 (2008) 2889–2904. [19] V.I. Terekhov, V.V. Terekhov, Heat transfer in a high vertical enclosure with fins attached to one of the side walls, High Temp. 44 (2006) 436–441. [20] L.A. Moiseeva, S.G. Cherkasov, Steady-state natural-convection heat transfer in a cylindrical tank under conditions of uniform heat supply and simultaneous heat removal through local sinks, High Temp. 35 (1997) 553–558. [21] G.V. Kuznetsov, M.A. Sheremet, The Rayleigh–Benard instability in an enclosure having finite thickness walls, J. Phys: Conf. Ser. 216 (2010) 1–15. [22] L. Valencia, J. Pallares, I. Cuesta, F. Xavier Grau, Turbulent Rayleigh–Benard convection of water in cubical cavities: a numerical and experimental study, Int. J. Heat Mass Transfer 50 (2007) 3203–3215. [23] M.Y. Ha, M.J. Jung, A numerical study on three-dimensional conjugate heat transfer of natural convection and conduction in a differentially heated cubic enclosure with a heat-generating cubic conducting body, Int. J. Heat Mass Transfer 43 (2000) 4229–4248. [24] Y.H. Li, K.C. Lin, T.F. Lin, Computation of unstable liquid metal convection in a vertical closed cylinder heated from the side and cooled from above, Numer. Heat Transfer Part A: Appl. 32 (1997) 289–309. [25] S.S. Leong, Numerical study of Rayleigh–Benard convection in a cylinder, Numer. Heat Transfer Part A: Appl. 41 (2002) 673–683. [26] T.C. Cheng, Y.H. Li, T.F. Lin, Effects of thermal boundary condition on buoyancy driven transitional air flow in a vertical cylinder heated from below, Numer. Heat Transfer Part A: Appl. 37 (2000) 917–936. [27] L.G. Loicyanski, Classical Fluid Mechanics, Drofa, Moscow, 2003 (p. 840). [28] G.V. Kuznetsov, M.A. Sheremet, Conjugate natural convection in an enclosure with local heat sources, Comput. Therm. Sci. 1 (2009) 341–360. [29] K. Aziz, J.D. Hellums, Numerical solution of three-dimensional equations of motion for laminar natural convection, Phys. Fluids 10 (1967) 314–324. [30] G.J. Hirasaki, J.D. Hellums, A general formulation of the boundary conditions on the vector potential in three dimensional hydrodynamics, Quart. Appl. Math. 16 (1968). [31] V.I. Paasonen, Boundary conditions of high-order accuracy at the poles of curvilinear coordinate systems, Russ. J. Numer. Anal. Math. Model. 14 (1999) 369–382. [32] A.A. Samarskii, Theory of Difference Schemes, Nauka, Moscow, 1977. 656 p. [33] V.M. Paskonov, V.I. Polezhaev, L.A. Chudov, Numerical Modelling of Heat and Mass Exchange Processes, Nauka, Moscow, 1984. 288 p. [34] M.A. Sheremet, Three-dimensional conjugate natural convection in a vertical cylinder under heat transfer to the surroundings, Fluid Dyn. 46 (2011) 647– 657. [35] B. Gebhart, Y. Jaluria, R.L. Mahajan, B. Sammakia, Buoyancy – Induced Flows and Transport, Hemisphere Publishing Corporation, New York, 1988. 1001 p.