- Email: [email protected]
A numerical study of three-dimensional laminar mixed convection past an open cavity
International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
A numerical study of three-dimensional laminar mixed convection past an open cavity Y. Stiriba *, F.X. Grau, J.A. Ferré, A. Vernet Universitat Rovira i Virgili, ETSEQ-DEM, Department of Mechanical Engineering, Av. Paisos Catalans 26, 43007 Tarragona, Spain
a r t i c l e
i n f o
Article history: Available online 1 July 2010 Keywords: Mixed convection Direct numerical simulation Open cavity Stability
a b s t r a c t A numerical study has been carried out to analyze the effects of mixed convective flow over a threedimensional cavity that lies at the bottom of a horizontal channel. The vertical walls of the cavity are isothermal and all other walls are adiabatic. The cavity is assumed to be cubic in geometry and the flow is laminar and incompressible. A direct numerical simulation is undertaken to investigate the flow structure, the heat transfer characteristics and the complex interaction between the induced stream flow at ambient temperature and the buoyancy-induced flow from the heated wall over a wide range of the Grashof number (103–106) and two Reynolds numbers Re = 100 and 1000. The computed thermal and flow fields are displayed and discussed in terms of the velocity fields, streamlines, the temperature distribution and the averaged Nusselt number at the heated and cooled walls. It is found that the flow becomes stable at moderate Grashof number and exhibit a three-dimensional structure, while for both high Reynolds and Grashof numbers the mixed convection effects come into play, push the recirculating zone further upstream and the flow becomes unsteady with Kelvin–Helmholtz instabilities at the shear layer. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Mixed convection and heat transfer over open cavities have been the subject of extensive research for many years due to their importance in various engineering systems. For instance, solar energy collectors, air cooling systems in electronic components, buildings, landing systems of aircrafts, nuclear reactors are just few examples. The presence of such type of geometry is of interest because of the convective heat transfer that occurs between the cavity and the forced flow stream of air (e.g., electronic devices, see [8]). Therefore a qualitative characterization of the complex interaction between the natural convection due to buoyancy forces and the external flow is highly important in the system design. Numerous experimental and numerical studies of flow past open cavity have been reported for both natural and mixed convections, see for instance [1–3,5,10,13,15] and references therein. In [10,12] the authors numerically studied the transport process occurring due to the interaction of air streams with buoyancy-induced flow within rectangular enclosures with constant source of heat for different flow parameters. Their results showed that for high value of the Richardson number (Ri) the average heat transfer rate increases at a fixed Reynolds number (Re). Beyond a critical value of the Richardson number the numerical solution gives unstable results. More detailed numerical studies of this oscillatory * Corresponding author. E-mail address: [email protected] (Y. Stiriba). 0017-9310/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2010.06.012
behavior was reported in [11]. A numerical investigation of the two-dimensional mixed heat transfer over vertical open cavity was carried out in [6] with three different angles of attack, heated horizontal walls, different Rayleigh numbers between 102 and 105, Reynolds numbers varying from 102 to 104 and aspect ratios varying from 0.25 to 1.0. The results indicated that the average Nusselt number increases on the upper and lower walls almost linearly, up to some critical Reynolds number, with the Reynolds and Grashof numbers and the horizontal flow can insulate the cavity from external flow. Manza et al. [7] investigated numerically, the mixed steady state convection in a two-dimensional horizontal open cavity using a Galerkin finite element method. They obtained results by uniformly heating the three walls separately for Richardson number from 0.1 to 100, Reynolds number 100–1000, and aspect ratio 0.1–1.5 and showed that the temperature increases with the Grashof number, and the aspect ratio has a significant effect on the Nusselt number, streamlines and isotherms. It was found that the Nusselt number and the maximum temperature were higher when the heated wall is on the opposite flow side. Zdanski et al. [19] presented numerical results for both laminar and turbulent flows over shallow cavities by varying the Re number, the cavity aspect ratio and the turbulence level of the incoming flow. They showed that the maximum heat transfer at the floor of the cavity coincides with the position where the turbulent diffusion position close to the wall is also a maximum, and the downstream step region is unaffected by the level of turbulence. Vafai and Huang [16] studied the effects of using intermittency porous open cavities for
4798
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
regulating and modifying the flow and heat transfer. By their configuration they analyzed the interaction phenomena occurring in the porous medium and the fluid layer. Recently Manza et al. [8] have investigated experimentally mixed convection in open cavities with heated inflow wall and Richardson numbers varying from 30 to 110 for Re = 1000 and 2800–8700 for Re = 100 and concluded that the surface temperature is lower for low Reynolds numbers at the some Richardson number. They observed also, for Re = 1000, a parallel forced flow motion and a recirculation flow inside the cavity. The two-dimensional version of this problem has received considerable attention. However, to the best of our knowledge, very few results have been obtained for the three-dimensional case. In this work we present a numerical study of thermal convection in a three-dimensional open cavity at Prandtl number Pr = 0.71. The cavity is located at the bottom of a horizontal channel. A parametric study is performed for Grashof number from 103 to 106 and two Reynolds numbers, Re = 100 and 1000, to investigate the flow structure and the heat transfer. The thermal and flow fields are presented and analyzed in terms of the velocity fields, streamlines and the temperature distribution. The interaction between the forced flow at ambient temperature and the buoyancy-induced flow from the heated wall is studied. Averaged Nusselt numbers at the heated and cooled walls for different flow parameters are also reported.
H
U0 L
g
T
y
C
TH D
x z
L
Fig. 1. Sketch of the computational domain.
2. Physical problem and governing equations
Fig. 2. Isotherms obtained for Re = 100, Gr = 104 (left) and Gr = 105 (right).
Table 1 Comparison of the Nusselt number obtained using the present code (Right) with those of Khanafer et al. [6]. Values of reference [6] were approximately read. Re = 1500
Nu (upper) Nu (lower)
Re = 2000
Ref. [6]
Present
Ref. [6]
Present
3.76 4.85
3.93 4.94
2.00 2.00
1.910 1.98
The overall geometry and the computational domain of the open cavity used in this study are shown in Fig. 1. The cavity has length L, width L, and depth D. Here we consider a cubic cavity (L/D = 1), which is similar to the bi-dimensional configuration used by Manza et al. [7]. The length of the bottom wall of the channel behind the cavity in the streamwise direction is L, and the height is H = L/2. The right-hand side wall of the cavity in the outflow direction is assumed to be at constant temperature TH and the left-hand side wall of the cavity is maintained at ambient temperature T1 < TH. The remaining walls are considered as adiabatic ones. The surrounding fluid interacting with the cavity is at temperature T1. The cooling flow is air having a Prandtl number Pr of 0.71 and is assumed to be laminar and incompressible. The
t
t
0.9
0.9
0.75
0.75
0.6
0.6
0.45
0.45
0.3
0.3
0.15
0.15
0
0
t
t
0.9
0.9
0.75
0.75
0.6
0.6
0.45
0.45
0.3
0.3
0.15
0.15
0
0
Fig. 3. Temperature distribution for Re = 100 at the z-central plan. (a) Ri = 0.1, (b) Ri = 1 (c) Ri = 10 and (d) Ri = 100.
4799
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
viscous dissipation is negligibly small and all other fluid properties are assumed constant except the fluid density in the buoyancy terms according to the Boussinesq approximation. Under the above assumptions, the continuity, momentum and energy equations for a three-dimensional laminar incompressible fluid can be expressed in non-dimensional form as follows:
@ui ¼0 @xi @ui @ðuj ui Þ @p 1 @ 2 ui Gr ¼ þ þ T di2 þ @t @xi Re @xj xj Re2 @xj
ð1Þ ð2Þ
@T @ðuj T Þ 1 1 @2T þ ¼ Re Pr @xj xj @t @xj
ð3Þ
with the following dimensionless variables
xi ¼
xi ; L
t ¼
U0t ; L
ui ¼
ui ; U0
p ¼
p p1 2 1 U0
q
;
T ¼
T T1 TH T1
ð4Þ
where p1 and q1 are reference pressure and density respectively, TH is the temperature at the heated surface, ui are the velocity com-
ponents, p is the pressure, Re = (U0L)/m is the Reynolds number, Gr = (gbDTL3)/m2 is the Grashof number, and Pr = m/a is the Prandtl number. Here b, m and a are the coefficient of volumetric expansion, kinematic viscosity, and thermal diffusivity, respectively. The extension of the open-ended domain requires substantially larger memory and computational time. Vafai and Ettefagh, see [13–15], checked the influence of outer boundary conditions on the flow and the temperature fields through various numerical runs covering different flow parameters such as the Raleigh number. They showed the validity of using symmetric conditions at the centerline of the cavity and far field conditions at other parts of the domain. Their approximations resulted in an accurate description of the flow inside the cavity. To save memory requirement and CPU time in three-dimensions, we used the same boundary conditions except at the bottom of the cavity. In relation to Fig. 1, the boundary conditions used are the following: at the inlet, the flow has uniform velocity U0 in the free stream direction:
u ¼ 1;
v ¼ 0;
w ¼ 0;
T ¼ 0;
x ¼ 0;
y P D=L
Fig. 4. Temperature distribution for Re = 100 at different x- and y-constants plans. Top: Ri = 0.1, center: Ri = 10 and bottom: Ri = 100.
4800
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
at the top of the channel, outflow condition was used:
@ui @y
¼ 0;
@T ¼ 0; @y
y ¼
DþH L
at the outlet, a boundary condition of convective type can be applied, which consists in solving the equation: @ui @ui ¼ 0; þ ui @t @xi
@T ¼ 0; @xi
x ¼ 2;
y P
D L
ð5Þ
at each time iteration, and the temperature has zero gradient. No-slip boundary condition is used at the remaining walls. The local Nusselt number (NuL) and surface averaged Nusselt number (NuS) are used to represent heat transfer at the heated and cooled walls
NuL ¼
@T ; @x x ¼0;1
NuS ¼
Z 0
L
Z
L
NuL dy dz :
ð6Þ
0
channel with vertical flow investigated by Khanafer et al. [6]. We use the same boundary conditions and wall temperatures on the horizontal walls of the cavity. Isotherm lines for Re = 100 and Gr = 104, 105 are presented in Fig. 2. It can be seen from Fig. 2 that the present results and those reported in [6] are in excellent agreement. The comparison of the Nusselt number for Re = 1500 and 2000 was shown in Table 1. The present solutions obtained by similar simulation fairly agree with the predicted results in reference [6]. A systematic grid independence study is carried out to check the grid independence of the solutions. Three different non-uniform grid sizes (120 90 60, 100 75 50, and 80 60 40), with a concentration of grid lines near the walls of the cavity, were used and finally (100 75 50) grid points were adopted for grid-free solution throughout the calculations in the present study. Note that the above tested grid sizes correspond to 603, 503 and 403 grid points inside the enclosure, respectively. The difference between the averaged surface Nusselt numbers corresponding to grid sizes 603 and 503 inside the cavity is 1.5%.
3. Numerical solution method
4. Results and discussion
In the present work, a parallel finite volume code is used. The flow solver discretizes the time-dependent incompressible Navier–Stokes Eqs. (1)–(3) in three-dimensions on non-uniform, staggered Cartesian meshes. The convective fluxes across the surfaces of the control volume are approximated using the sharp and monotonic algorithm for realistic transport (SMART) scheme [4], and the diffusive terms with central differences. The numerical scheme employs the Adams–Bashforth method for time marching. The continuity and momentum equations are coupled via the SMACmethod, in which, the Poisson equation for the pressure is computed with the biconjugate gradient method (BiCGtab) [17]. The code is parallelized using the domain decomposition method and Message Interface Passing (MPI) libraries to communicate data between processors. The code was validated for several problems including lid driven cavity and natural convection in cubical cavities [9]. Here we concentrate only on those tests which involve mixed convection in open closures, namely the two-dimensional open cavity in a
In order to study the mixed convection flow effects, the computed thermal and flow fields are analyzed in terms of the temperature distribution, the development of the velocity components and the heat transfer at the heated wall for different Ri = Gr/Re2 values between 103 and 102. The Reynolds numbers examined were 100 and 1000, for which the flow is expected to remain in the laminar regime. Papanicolaou and Jaluria [10] observed that increasing further the Re number the solution did not converge and a turbulent closure model is required to simulate the flow. 4.1. Temperature profile Computed temperature distributions are presented along different plans in Figs. 3–6 for different Ri numbers. Fig. 3 shows the temperature contours at the z-central plan obtained by varying the Ri number, or equivalently Gr, at Re = 100. As discussed earlier, see [7], for these ranges the steady state situation is maintained. For low Ri number values, see Figs. 3(a) and (b) and 4, the
Fig. 5. Temperature distribution for Re = 1000 at the z-central plan. (a) Ri = 0.01, (b) Ri = 0.1 (c) Ri = 1 and (d) Ri = 10.
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
4801
Fig. 6. Temperature distribution for Re = 1000 at different x- and y-constants plans. Top: Ri = 0.01, center: Ri = 0.1 and bottom: Ri = 1.
isothermlines are clustered close to the top right wall and are parallel to the heated and cooled walls. The entire region of the cavity is heated and the cold external air flow passes through the cavity without penetrating downward. However, the heated zone of the enclosure becomes smaller as the Ri number increases to 10 and 100, and the isotherms are more and more concentrated at the lower right corner. This can also be observed in Fig. 4 where the isotherms in the constants x- and y-plans form a layer along the heated wall. This means that high temperature gradients are present near the heated and cooled walls due to the imposed boundary conditions, and the heat is transferred by conduction for low Ri values. For higher Richardson number of Ri = 100, the isothermlines are almost parallel to the right vertical wall and uniformly narrow. Furthermore, the heating becomes stronger along the bottom horizontal wall of the channel in the streamwise direction and the rest of the cavity is almost unaffected, see also Fig. 4. The outgoing fluid rises up at much faster speed. This indicates that the buoyancy effects are much stronger than that of the forced flow and forces the heated air to rise up in the direction which is opposite to that of the gravitational field.
This fluid mechanics behavior is in agreement with previous results obtained using two-dimensional physical model by Manza et al. [7]. Nevertheless, the flow has a three-dimensional nature as we can observe from Fig. 4. The same ‘ejection mechanism’ were observed in [13,15] with a thermal boundary layer along the heated walls and outside the enclosure. Fig. 5 shows the temperature profiles at the z-central plan for Re = 1000. For Ri = 0.01 and Ri = 0.1 the solution reaches a steady state. For low Ri values, the isothermlines are more clustered at the top right corner vertical wall and parallel to the heated and cooled walls. The heat is transferred to the entire space of the cavity mainly by conduction. The upper part of the enclosure, right below the external flow, has horizontal stratified isotherms which means that the forced flow is stronger and the flow inside the cavity is dominated by heat conduction mechanism. As Ri increases to 0.1 the contour lines deviate from the middle heated wall to the upper part of the cavity due to the dominance of the natural convection effects. In this case the buoyancy forces become stronger and push the heated air further upstream. From Fig. 6 we observe
4802
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
Fig. 7. Flow structures for Re = 100 (top and middle) and Re = 1000 (bottom), (a) Ri = 0.1, (b) Ri = 1, (c) Ri = 10, (d) Ri = 100, (e) Ri = 0.01 and (f) Ri = 0.1.
that the flow exhibits a three-dimensional structure. As the Ri increases further (Ri P 1) the flow becomes unsteady. The temperature profile indicates the Kelvin–Helmholtz instabilities in the shear layer, at Ri = 1, due to the interaction between the forced flow and the buoyancy-induced flow by the natural convection. This instability is located right below the external flow close to the heated wall similar to the one noticed by Yao et al. [18] at high Re number without heat transfer. The flow becomes extremely unstable with unsteady solution. This is why the results are not presented here. Note that the isothermal wall remain parallel in the neighborhood of the cooled wall. At Ri = 10 the isothermlines are very close, the thermal boundary layer from the heated wall becomes thinner, and the hot buoyant fluid rises like a buoyant plume into the outside domain and escape from the restrictive vertical cavity walls. The thickness of this thermal layer decreases as Ri increases, as seen in Fig. 5(c) and (d) as in [13,15]. We found four grid nodes inside this thermal layer. The isotherm pattern indicates that up to a critical value of the Ri number, the buoyancy forces (temperature gradients) exceed the dissipative effects of viscous forces causing wavy pattern of the isotherms at the outflow channel. Its effect is to enhance heat transfer rate of the external flow leaving the channel and generate more unstable wakes. A thermal plume can also be seen to rise into the floor of the cavity, see Fig. 6.
4.2. Development of velocity components To see the impact of the mixed convection on the flow pattern and to give the idea of recirculating and ascending flow structure, the streamlines, and v, w-velocity profiles at specific planes are displayed in Figs. 7–11. In Fig. 7 we present the streamlines map for Re = 100. For Ri = 0.1 in Fig. 7(a), the cold forced flow passes through the cavity without penetrating and creates a recirculating structure that fills the entire cavity. As the Ri was increased, in addition to the primary eddy two other minor eddies can be seen at the left and right bottom corner, see Fig. 7(b), since the heated and cooled airs move upward and downward, respectively, due to the buoyancy forces. For Ri = 10 and 100, a recirculating cell is present at the upper left of the cavity. The temperature gradients becomes stronger and part of the heated flow leaves close to the heated wall and the horizontal adiabatic wall of the channel behind the enclosure. Part of the forced flow also penetrates further downstream in the enclosure. These explain the higher heated transfer already observed in the temperature profiles. Fig. 8 shows the transverse velocity component (v/U0) as a function of the streamwise distance (x/L) at different constant spanwise and vertical plans. Near the side wall (x/L = 0), and for low Ri numbers (see Fig. 8(a) and (b)), the v-velocity component profile presents positive values meaning that the fluid is moving up and
4803
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808 0.15
0.15
z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
0.1
0.1 0.05
v/U
0.05
v/U
z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
0
0
-0.05
-0.05
-0.1
-0.1
-0.15
-0.15 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
x/L 0.3
0.6
0.8
1
x/L z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
1.5
0.2 1
v/U
v/U
0.1
0
0.5
0
-0.1
-0.5
-0.2
-1 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
x/L
0.6
0.8
1
x/L
Fig. 8. v-velocity component for Re = 100 at x- and y-planes, (a) Ri = 0.1, (b) Ri = 1, (c) Ri = 10 and (d) Ri = 100.
0.06
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.04 0.02
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.04 0.02 0
-0.02
v/U
v/U
0
-0.02
-0.04
-0.04
-0.06
-0.06
-0.08
-0.08
-0.1
-0.1 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
z/L 1
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.3
0.6
0.8
1
z/L x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.8 0.6
0.2
v/U
v/U
0.4 0.1
0.2 0
0
-0.2 -0.1 -0.4 -0.2
-0.6 0
0.2
0.4
0.6
z/L
0.8
1
0
0.2
0.4
0.6
z/L
Fig. 9. v-velocity component for Re = 100 at y- and z-planes, (a) Ri = 0.1, (b) Ri = 1, (c) Ri = 10 and (d) Ri = 100.
0.8
1
4804
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
0.06
0.06
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.04
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.04
0.02
w/U
w/U
0.02
0
0
-0.02
-0.02
-0.04
-0.04 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
z/L 0.2
0.3
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.15
0.6
0.8
1
z/L
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.2
0.1
0.05
w/U
w/U
0.1
0
0
-0.05 -0.1 -0.1 -0.15
-0.2 0
0.2
0.4
0.6
0.8
1
z/L
0
0.2
0.4
0.6
0.8
1
z/L
Fig. 10. w-velocity component for Re = 100 at y- and z-planes, (a) Ri = 0.1, (b) Ri = 1, (c) Ri = 10 and (d) Ri = 100.
negative values near the side wall (x/L = 1) meaning that the flow is moving down up. This leads to the primary recirculation, similar to the natural convection profiles, already observed in the recirculating nature in the flow structure pictures. Near the bottom wall (y/L = 0.25) the convective flow remains almost constant for low Ri values, see Fig. 8(a) and (b). Fig. 8(b) indicates a slight increase in the v-component, for Ri = 0.1, near the side wall (x/L = 1) meaning a quasi-two-dimensional secondary eddy. At Ri = 1 and 10, see Fig. 8(c) and (d), the buoyancy effects considerably modify the v-velocity. The variation is much steeper near the heated wall (x/L = 1). The heated air rises fast with increasing Ri. For Ri = 1, the change of the v-velocity is small at the bottom and middle part near the side wall (x/L = 0), as opposed to the top left side of the cavity, where the v-component experiences a change in the slop. At the bottom left side of the enclosure the v-component presents negative values meaning that the flow is sinking along the cooled wall. At Ri = 10, the vertical velocity is almost constant near the side wall (x/L = 0). The distribution gives the idea of recirculating and ascending flow structure. Similar behaviors were observed in the two-dimensional results of Manza et al. [7] at Ri = 0.1 and Ri = 100 for Re = 100. Fig. 9 shows the v-velocity component as a function of the spanwise distance (z/L) at two different planes, the constant plan x/L = 0.5 at the middle of the cavity and the plan x/L = 0.75 close to the step. The value y/L = 0.75 corresponds to a transverse plan near the top enclosure, while y/L = 0.5 corresponds to a transverse plan at the middle cavity. At the middle constant (x/L = 0.5) plan and for both Ri = 0.01 and 0.1, the v-velocity distribution shows that the flow remains nearly constant. It descends toward the bottom as Ri increases to 1 and 10. Close to the step, the flow ascends toward the top in the central direction for Ri = 0.1 and descends
downward for high Ri values. Near the top enclosure and the heated wall the flow is displaced towards the bottom in the central direction for Ri = 0.01 and 0.1. At Ri = 1 and 10, the vertical velocity component behavior is similar to Fig. 9(a) and (b) discussed above, the flow is ascending in the vicinity of the side walls meaning that the hot air is rising, see also Fig. 6. Fig. 10 shows the w-velocity as a function of the spanwise distance (z/L) at the same constant plans as in Fig. 9. At the middle plan, the flow is directed toward the center of the enclosure. Near the step, the same figure shows that the flow is moving toward the side walls. Even at low Ri, where the flow is quasi-two-dimensional, the movement is always three-dimensional with a non negligible w-component (see Fig. 10(a) and (b)). Flow structures for Ri = 0.01 and 0.1 at Re = 1000 are depicted in Fig. 7(c) and (d). For Ri = 0.01 the streamlines have the same behavior as for Re = 100. The amount of recirculation flow increases due to stronger forced flow effects. The vertical velocity profile in Fig. 11(c) highlights the recirculating cell. For Ri = 0.1 one can observe the development of two recirculating structures. One cell is located just below the airstream and the other cell is in the downstream of the enclosure due to the raised heated air and the sinking cooled air. Both eddies are accompanied with a stratified air at the floor of the cavity just below streamlines of external flow. The velocity component profiles in Fig. 11(d) verify this scenario. Fig. 11 shows also the v- and w-velocity components as a function of the spanwise distance (z/L) and streamwise distance (x/L) at two at two different planes, at the middle and close to the floor cavity. It can be seen that the vertical velocity v, see Fig. 11(a) and (b), increases in magnitude because the forced flow moves faster with increasing Re. The velocities inside the enclosure experiences irregular ascending and descending at the vertical and
4805
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
0.15
0.15
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.1
0.1 0.05
0
v/U
v/U
0.05
-0.05
0 -0.05
-0.1
-0.1
-0.15 -0.2
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0
0.2
0.4
0.6
0.8
-0.15
1
0
0.2
0.4
z/L 0.15
0.2
z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
0.1
1
0.1
0 -0.05
v/U
v/U
0.8
z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
0.15
0.05
-0.1 -0.15
0.05 0 -0.05
-0.2
-0.1
-0.25 -0.3
0.6
z/L
0
0.2
0.4
0.6
0.8
-0.15
1
0
0.2
0.4
x/L 0.06
0.08
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.04
0.6
0.8
1
x/L
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.06 0.04
w/U
w/U
0.02
0
0.02 0 -0.02
-0.02
-0.04 -0.04
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
z/L
0.6
0.8
1
z/L
Fig. 11. Velocity components for Re = 1000, left: Ri = 0.01 and right: Ri = 0.1.
horizontal directions and the fluid is no longer steady, see also Fig. 11(f). The w-velocity component, Fig. 11(e), displays the same behavior as in the previous case with Re = 100. The effect of varying the ratio of the height of the vertical heated wall D is shown in Figs. 12(a) and (c) and 13(a) and (c) for Re = 100 and Ri = 0.1. Figs. 12(a) and (c) and 13(a) and (c) show the isotherms and streamlines, respectively, for three cavities of aspect ratio AR = 0.5, 1.5 and 2. For AR = 0.5, the main vortex fills the cavity close to the opening vertical wall. For a deeper cavity with AR = 1.5, the isotherms remains parallel to the heated wall due to the high gradients in the vertical direction. In addition, the intensity of the cavity flow decreases at the bottom wall and the fluid flows along the lower plate. For AR = 2, a pair of counter-rotating vortices appears in the cavity. This can be attributed to the buoyancy forces which opposes the recirculating flow within
the enclosure. The fluid in this case penetrates into the cavity and reverses its direction in the bottom region. The isotherms in the lower part are almost parallel to the vertical heated wall. The isotherms and streamlines for high Richardson number of Ri = 1 and aspect ratio AR = 2 are presented in Figs. 12(d) and 13(d). As the Ri increases, the second vortex which is caused by the strong buoyancy forces becomes much stronger and occupies more than the half of the open cavity. This recirculating flow is responsible for the penetration and inversion of the fluid at the central half portion of the cavity. The time behavior of flow variables at the trailing edge corner for Ri = 1 is shown in Fig. 14. The flow enters the transient regime and becomes unstable. For both high Re and Ri numbers the buoyancy forces exceed the dissipative effects of viscous forces causing the fluid to rise quickly and forming an unsteady convection and fluctu-
4806
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
number at the hot wall is always seen to increase with Ri for each value of Re. Similar behavior were observed in [10] by the author in their bi-dimensional numerical tests. At low values of Ri number and for both Re = 100 and Re = 1000, the Nusselt number increases slightly since the buoyancy effects are weak and the mixed convection transport mechanism is largely dominated by convective effects. Significant increase is found for high values of Ri, up to certain Ri value for each Re number, since the natural convection was the dominant mode of the mixed convection transport. The continuous increase of the heat transfer with increasing Reynolds for each Ri value is also observed. 5. Conclusions
Fig. 12. Temperature distribution for Re = 100 and Ri = 0.1, (a) D = 0.5, (b) D = 1.5, (c) D = 2, and (d) D = 2 with Ri = 1.
ation in velocity and temperature. This behavior will need more detailed study and would constitute a subject of a future work. 4.3. Heat transfer The averaged Nusselt number calculated at the hot and cooled wall is presented against Ri for Re = 100 and 1000 in Fig. 15. At cooled wall, the Nusslet number decreases slightly with Ri. It becomes very small for Re = 100 and Ri = 100 since the isotherm lines are clustered parallel to the hot wall, see also Fig. 4. The Nusselt
A detailed numerical study has been carried out to investigate mixed convection for incompressible laminar flow past open cubical cavity. The results show that the flow exhibits a three-dimensional structure and becomes steady for Re = 100 with Ri 6 100 and Re = 1000 with Ri 6 0.1. The forced flow dominates the flow transport mechanism and large recirculating zone form inside the enclosure which results in heat transfer by conduction. In this case the Nusselt number increases slightly. For both high Reynolds and Richardson numbers the natural convection comes into play and push the recirculating zone further upstream. The flow becomes unsteady with Kelvin–Helmholtz instabilities at the shear layer and the Nusselt number increases significantly. For Re = 1000 and Gr = 107 the isotherm lines remain parallel close to the heated wall and form and thermal boundary layer. The work also gave quantitative results for the velocity components resulting from various locations to give the idea of recirculating and ascending flow structure. This flow has revealed to be a very convenient case for studying flow transients under the combined effects of forced and natural convection with a progressive evolution from a quasi-two-dimensional flow to a three-dimensional steady and
Fig. 13. Flow structures for Re = 100 and Ri = 0.1, (a) D = 0.5, (b) D = 1.5, (c) D = 2 and (d) D = 2 with Ri = 1.
4807
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
0.6
0.6 0.5
0.5
0.4 0.4 0.3 0.2
V
U
0.3 0.2
0.1 0
0.1
-0.1 0 -0.2 -0.1 -0.2
-0.3 20
40
60
80
100
120
140
160
-0.4
20
40
60
Time
80
100
120
140
160
Time
0.3
0.5
0.2
0.4
0.1
T
W
0.3 0
0.2 -0.1
0.1
-0.2
-0.3
20
40
60
80
100
120
140
160
0
20
40
Time
60
80
100
120
140
160
Time
Fig. 14. Time history of the instantaneous velocity components and temperature for Re = 1000 and Ri = 1.
16
Acknowledgments
14
This work was financially supported by project DPI2006-02477 from Ministerio de Educación y Ciencia (Spain) and Fondos FEDER.
12 10
References
Nu
8 6 4 2 0 -2 -4 10
-4
10
-3
10
-2
10
-1
10
0
10
1
10
2
10
3
Ri Fig. 15. Average Nusselt number vs. Ri numbers. Cold wall, ‘D’: Re = 100 and ‘r’: Re = 1000. Hot wall, ‘’: Re = 100 and ‘}’: Re = 1000.
unsteady flow in a relatively tidy range of Re and Ri numbers. Future work are being undertaken to study the unsteady behavior for both high Re and Ri.
[1] E. Bilgen, H. Oztop, Natural convection heat transfer in partially open inclined square cavities, Int. J. Heat Mass Transfer 48 (2005) 1470–1479. [2] E. Chakroun, M.M. Elsayed, S.F. Al-Fahed, Experimental measurements of heat transfer coefficient in partially/fully opened titled partially open cavities, J. Solar Eng. 119 (1997) 298–303. [3] Y.L. Chan, C.L. Tien, A numerical study of two dimensional natural convection in square cavities, Numer. Heat Transfer 8 (1985) 65–80. [4] P. Gaskell, A. Lau, Curvature-compensated convective transport: SMART, a new boundedness-preserving transport algorithm, Int. J. Numer. Meth. Fluids 8 (1988) 617–641. [5] A. Javam, S.W. Armfield, Stability and transition of stratified natural convection flow in open cavities, J. Fluid Mech. 445 (2001) 285–303. [6] K. Khanafer, K. Vafai, M. Lightstone, Mixed convection heat transfer in twodimensional open-ended enclosures, Int. J. Heat Mass Transfer 45 (2002) 5171–5190. [7] O. Manza, S. Nardini, K. Khanafer, K. Vafai, Effect of heated wall position on mixed convection with an open cavity, Numer. Heat Transfer, Part A 43 (2003) 259–282. [8] O. Manza, S. Nardini, K. Vafai, Experimental investigation of mixed convection in a channel with an open cavity, Exp. Heat Transfer 19 (2006) 53–68. [9] J. Pallares, I. Cuesta, F.X. Grau, F. Giralt, Natural convection in a cubical cavity heated from bellow at low Rayleigh numbers, Int. J. Heat Mass Transfer 39 (1996) 3233–3247.
4808
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
[10] E. Papanicolaou, Y. Jaluria, Mixed convection from an isolated heat source in a rectangular enclosure, Numer. Heat Transfer, Part A: Appl. 18 (4) (1990) 427– 461. [11] E. Papanicolaou, Y. Jaluria, Transition to a periodic regime in mixed convection in a square cavity, J. Fluid Mech. 239 (1992) 489–509. [12] E. Papanicolaou, Y. Jaluria, Mixed convection from an localized heat source in a cavity with conducting walls, Numer. Heat Transfer, Part A: Applications 23 (4) (1990) 463–484. [13] K. Vafai, J. Ettefagh, The effects of sharp corners on buoyancy-driven flow with particular emphasis on outer boundaries, Int. J. Heat Mass Transfer 33 (1990) 2318–2328. [14] K. Vafai, J. Ettefagh, Thermal and fluid flow instabilities in buoyancy-driven flows in open-ended cavities, Int. J. Heat Mass Transfer 33 (1990) 2329–2344.
[15] K. Vafai, J. Ettefagh, Thermal and fluid flow instabilities in buoyancy-driven flows in open-ended cavities, ASME J. Heat Transfer 116 (1994) 604–613. [16] K. Vafai, P.C. Huang, Analysis of heat transfer regulation and modification employing intermittency emplaced porous cavities, Int. J. Heat Mass Transfer 33 (1990) 2329–2344. [17] H.A. Van der Vost, BiCGSTAB: a fast and smooth converging variant of Bi-CG for the solution of non-symmetric linear system, SIAM J. Sci. Statist. Comput. 13 (1992) 631–644. [18] H. Yao, R.K. Cooper, S. Raghunathan, Numerical simulation of incompressible laminar flow over three-dimensional rectangular cavities, J. Fluid Eng. 126 (2004) 919. [19] P.S.B. Zdanski, M.A. Ortega Jr., N.G.C.R. Fico, Heat transfer in the flow over shallow cavities, J. Heat Transfer ASME 127 (2005) 699–712.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
A numerical study of three-dimensional laminar mixed convection past an open cavity Y. Stiriba *, F.X. Grau, J.A. Ferré, A. Vernet Universitat Rovira i Virgili, ETSEQ-DEM, Department of Mechanical Engineering, Av. Paisos Catalans 26, 43007 Tarragona, Spain
a r t i c l e
i n f o
Article history: Available online 1 July 2010 Keywords: Mixed convection Direct numerical simulation Open cavity Stability
a b s t r a c t A numerical study has been carried out to analyze the effects of mixed convective flow over a threedimensional cavity that lies at the bottom of a horizontal channel. The vertical walls of the cavity are isothermal and all other walls are adiabatic. The cavity is assumed to be cubic in geometry and the flow is laminar and incompressible. A direct numerical simulation is undertaken to investigate the flow structure, the heat transfer characteristics and the complex interaction between the induced stream flow at ambient temperature and the buoyancy-induced flow from the heated wall over a wide range of the Grashof number (103–106) and two Reynolds numbers Re = 100 and 1000. The computed thermal and flow fields are displayed and discussed in terms of the velocity fields, streamlines, the temperature distribution and the averaged Nusselt number at the heated and cooled walls. It is found that the flow becomes stable at moderate Grashof number and exhibit a three-dimensional structure, while for both high Reynolds and Grashof numbers the mixed convection effects come into play, push the recirculating zone further upstream and the flow becomes unsteady with Kelvin–Helmholtz instabilities at the shear layer. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Mixed convection and heat transfer over open cavities have been the subject of extensive research for many years due to their importance in various engineering systems. For instance, solar energy collectors, air cooling systems in electronic components, buildings, landing systems of aircrafts, nuclear reactors are just few examples. The presence of such type of geometry is of interest because of the convective heat transfer that occurs between the cavity and the forced flow stream of air (e.g., electronic devices, see [8]). Therefore a qualitative characterization of the complex interaction between the natural convection due to buoyancy forces and the external flow is highly important in the system design. Numerous experimental and numerical studies of flow past open cavity have been reported for both natural and mixed convections, see for instance [1–3,5,10,13,15] and references therein. In [10,12] the authors numerically studied the transport process occurring due to the interaction of air streams with buoyancy-induced flow within rectangular enclosures with constant source of heat for different flow parameters. Their results showed that for high value of the Richardson number (Ri) the average heat transfer rate increases at a fixed Reynolds number (Re). Beyond a critical value of the Richardson number the numerical solution gives unstable results. More detailed numerical studies of this oscillatory * Corresponding author. E-mail address: [email protected] (Y. Stiriba). 0017-9310/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2010.06.012
behavior was reported in [11]. A numerical investigation of the two-dimensional mixed heat transfer over vertical open cavity was carried out in [6] with three different angles of attack, heated horizontal walls, different Rayleigh numbers between 102 and 105, Reynolds numbers varying from 102 to 104 and aspect ratios varying from 0.25 to 1.0. The results indicated that the average Nusselt number increases on the upper and lower walls almost linearly, up to some critical Reynolds number, with the Reynolds and Grashof numbers and the horizontal flow can insulate the cavity from external flow. Manza et al. [7] investigated numerically, the mixed steady state convection in a two-dimensional horizontal open cavity using a Galerkin finite element method. They obtained results by uniformly heating the three walls separately for Richardson number from 0.1 to 100, Reynolds number 100–1000, and aspect ratio 0.1–1.5 and showed that the temperature increases with the Grashof number, and the aspect ratio has a significant effect on the Nusselt number, streamlines and isotherms. It was found that the Nusselt number and the maximum temperature were higher when the heated wall is on the opposite flow side. Zdanski et al. [19] presented numerical results for both laminar and turbulent flows over shallow cavities by varying the Re number, the cavity aspect ratio and the turbulence level of the incoming flow. They showed that the maximum heat transfer at the floor of the cavity coincides with the position where the turbulent diffusion position close to the wall is also a maximum, and the downstream step region is unaffected by the level of turbulence. Vafai and Huang [16] studied the effects of using intermittency porous open cavities for
4798
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
regulating and modifying the flow and heat transfer. By their configuration they analyzed the interaction phenomena occurring in the porous medium and the fluid layer. Recently Manza et al. [8] have investigated experimentally mixed convection in open cavities with heated inflow wall and Richardson numbers varying from 30 to 110 for Re = 1000 and 2800–8700 for Re = 100 and concluded that the surface temperature is lower for low Reynolds numbers at the some Richardson number. They observed also, for Re = 1000, a parallel forced flow motion and a recirculation flow inside the cavity. The two-dimensional version of this problem has received considerable attention. However, to the best of our knowledge, very few results have been obtained for the three-dimensional case. In this work we present a numerical study of thermal convection in a three-dimensional open cavity at Prandtl number Pr = 0.71. The cavity is located at the bottom of a horizontal channel. A parametric study is performed for Grashof number from 103 to 106 and two Reynolds numbers, Re = 100 and 1000, to investigate the flow structure and the heat transfer. The thermal and flow fields are presented and analyzed in terms of the velocity fields, streamlines and the temperature distribution. The interaction between the forced flow at ambient temperature and the buoyancy-induced flow from the heated wall is studied. Averaged Nusselt numbers at the heated and cooled walls for different flow parameters are also reported.
H
U0 L
g
T
y
C
TH D
x z
L
Fig. 1. Sketch of the computational domain.
2. Physical problem and governing equations
Fig. 2. Isotherms obtained for Re = 100, Gr = 104 (left) and Gr = 105 (right).
Table 1 Comparison of the Nusselt number obtained using the present code (Right) with those of Khanafer et al. [6]. Values of reference [6] were approximately read. Re = 1500
Nu (upper) Nu (lower)
Re = 2000
Ref. [6]
Present
Ref. [6]
Present
3.76 4.85
3.93 4.94
2.00 2.00
1.910 1.98
The overall geometry and the computational domain of the open cavity used in this study are shown in Fig. 1. The cavity has length L, width L, and depth D. Here we consider a cubic cavity (L/D = 1), which is similar to the bi-dimensional configuration used by Manza et al. [7]. The length of the bottom wall of the channel behind the cavity in the streamwise direction is L, and the height is H = L/2. The right-hand side wall of the cavity in the outflow direction is assumed to be at constant temperature TH and the left-hand side wall of the cavity is maintained at ambient temperature T1 < TH. The remaining walls are considered as adiabatic ones. The surrounding fluid interacting with the cavity is at temperature T1. The cooling flow is air having a Prandtl number Pr of 0.71 and is assumed to be laminar and incompressible. The
t
t
0.9
0.9
0.75
0.75
0.6
0.6
0.45
0.45
0.3
0.3
0.15
0.15
0
0
t
t
0.9
0.9
0.75
0.75
0.6
0.6
0.45
0.45
0.3
0.3
0.15
0.15
0
0
Fig. 3. Temperature distribution for Re = 100 at the z-central plan. (a) Ri = 0.1, (b) Ri = 1 (c) Ri = 10 and (d) Ri = 100.
4799
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
viscous dissipation is negligibly small and all other fluid properties are assumed constant except the fluid density in the buoyancy terms according to the Boussinesq approximation. Under the above assumptions, the continuity, momentum and energy equations for a three-dimensional laminar incompressible fluid can be expressed in non-dimensional form as follows:
@ui ¼0 @xi @ui @ðuj ui Þ @p 1 @ 2 ui Gr ¼ þ þ T di2 þ @t @xi Re @xj xj Re2 @xj
ð1Þ ð2Þ
@T @ðuj T Þ 1 1 @2T þ ¼ Re Pr @xj xj @t @xj
ð3Þ
with the following dimensionless variables
xi ¼
xi ; L
t ¼
U0t ; L
ui ¼
ui ; U0
p ¼
p p1 2 1 U0
q
;
T ¼
T T1 TH T1
ð4Þ
where p1 and q1 are reference pressure and density respectively, TH is the temperature at the heated surface, ui are the velocity com-
ponents, p is the pressure, Re = (U0L)/m is the Reynolds number, Gr = (gbDTL3)/m2 is the Grashof number, and Pr = m/a is the Prandtl number. Here b, m and a are the coefficient of volumetric expansion, kinematic viscosity, and thermal diffusivity, respectively. The extension of the open-ended domain requires substantially larger memory and computational time. Vafai and Ettefagh, see [13–15], checked the influence of outer boundary conditions on the flow and the temperature fields through various numerical runs covering different flow parameters such as the Raleigh number. They showed the validity of using symmetric conditions at the centerline of the cavity and far field conditions at other parts of the domain. Their approximations resulted in an accurate description of the flow inside the cavity. To save memory requirement and CPU time in three-dimensions, we used the same boundary conditions except at the bottom of the cavity. In relation to Fig. 1, the boundary conditions used are the following: at the inlet, the flow has uniform velocity U0 in the free stream direction:
u ¼ 1;
v ¼ 0;
w ¼ 0;
T ¼ 0;
x ¼ 0;
y P D=L
Fig. 4. Temperature distribution for Re = 100 at different x- and y-constants plans. Top: Ri = 0.1, center: Ri = 10 and bottom: Ri = 100.
4800
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
at the top of the channel, outflow condition was used:
@ui @y
¼ 0;
@T ¼ 0; @y
y ¼
DþH L
at the outlet, a boundary condition of convective type can be applied, which consists in solving the equation: @ui @ui ¼ 0; þ ui @t @xi
@T ¼ 0; @xi
x ¼ 2;
y P
D L
ð5Þ
at each time iteration, and the temperature has zero gradient. No-slip boundary condition is used at the remaining walls. The local Nusselt number (NuL) and surface averaged Nusselt number (NuS) are used to represent heat transfer at the heated and cooled walls
NuL ¼
@T ; @x x ¼0;1
NuS ¼
Z 0
L
Z
L
NuL dy dz :
ð6Þ
0
channel with vertical flow investigated by Khanafer et al. [6]. We use the same boundary conditions and wall temperatures on the horizontal walls of the cavity. Isotherm lines for Re = 100 and Gr = 104, 105 are presented in Fig. 2. It can be seen from Fig. 2 that the present results and those reported in [6] are in excellent agreement. The comparison of the Nusselt number for Re = 1500 and 2000 was shown in Table 1. The present solutions obtained by similar simulation fairly agree with the predicted results in reference [6]. A systematic grid independence study is carried out to check the grid independence of the solutions. Three different non-uniform grid sizes (120 90 60, 100 75 50, and 80 60 40), with a concentration of grid lines near the walls of the cavity, were used and finally (100 75 50) grid points were adopted for grid-free solution throughout the calculations in the present study. Note that the above tested grid sizes correspond to 603, 503 and 403 grid points inside the enclosure, respectively. The difference between the averaged surface Nusselt numbers corresponding to grid sizes 603 and 503 inside the cavity is 1.5%.
3. Numerical solution method
4. Results and discussion
In the present work, a parallel finite volume code is used. The flow solver discretizes the time-dependent incompressible Navier–Stokes Eqs. (1)–(3) in three-dimensions on non-uniform, staggered Cartesian meshes. The convective fluxes across the surfaces of the control volume are approximated using the sharp and monotonic algorithm for realistic transport (SMART) scheme [4], and the diffusive terms with central differences. The numerical scheme employs the Adams–Bashforth method for time marching. The continuity and momentum equations are coupled via the SMACmethod, in which, the Poisson equation for the pressure is computed with the biconjugate gradient method (BiCGtab) [17]. The code is parallelized using the domain decomposition method and Message Interface Passing (MPI) libraries to communicate data between processors. The code was validated for several problems including lid driven cavity and natural convection in cubical cavities [9]. Here we concentrate only on those tests which involve mixed convection in open closures, namely the two-dimensional open cavity in a
In order to study the mixed convection flow effects, the computed thermal and flow fields are analyzed in terms of the temperature distribution, the development of the velocity components and the heat transfer at the heated wall for different Ri = Gr/Re2 values between 103 and 102. The Reynolds numbers examined were 100 and 1000, for which the flow is expected to remain in the laminar regime. Papanicolaou and Jaluria [10] observed that increasing further the Re number the solution did not converge and a turbulent closure model is required to simulate the flow. 4.1. Temperature profile Computed temperature distributions are presented along different plans in Figs. 3–6 for different Ri numbers. Fig. 3 shows the temperature contours at the z-central plan obtained by varying the Ri number, or equivalently Gr, at Re = 100. As discussed earlier, see [7], for these ranges the steady state situation is maintained. For low Ri number values, see Figs. 3(a) and (b) and 4, the
Fig. 5. Temperature distribution for Re = 1000 at the z-central plan. (a) Ri = 0.01, (b) Ri = 0.1 (c) Ri = 1 and (d) Ri = 10.
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
4801
Fig. 6. Temperature distribution for Re = 1000 at different x- and y-constants plans. Top: Ri = 0.01, center: Ri = 0.1 and bottom: Ri = 1.
isothermlines are clustered close to the top right wall and are parallel to the heated and cooled walls. The entire region of the cavity is heated and the cold external air flow passes through the cavity without penetrating downward. However, the heated zone of the enclosure becomes smaller as the Ri number increases to 10 and 100, and the isotherms are more and more concentrated at the lower right corner. This can also be observed in Fig. 4 where the isotherms in the constants x- and y-plans form a layer along the heated wall. This means that high temperature gradients are present near the heated and cooled walls due to the imposed boundary conditions, and the heat is transferred by conduction for low Ri values. For higher Richardson number of Ri = 100, the isothermlines are almost parallel to the right vertical wall and uniformly narrow. Furthermore, the heating becomes stronger along the bottom horizontal wall of the channel in the streamwise direction and the rest of the cavity is almost unaffected, see also Fig. 4. The outgoing fluid rises up at much faster speed. This indicates that the buoyancy effects are much stronger than that of the forced flow and forces the heated air to rise up in the direction which is opposite to that of the gravitational field.
This fluid mechanics behavior is in agreement with previous results obtained using two-dimensional physical model by Manza et al. [7]. Nevertheless, the flow has a three-dimensional nature as we can observe from Fig. 4. The same ‘ejection mechanism’ were observed in [13,15] with a thermal boundary layer along the heated walls and outside the enclosure. Fig. 5 shows the temperature profiles at the z-central plan for Re = 1000. For Ri = 0.01 and Ri = 0.1 the solution reaches a steady state. For low Ri values, the isothermlines are more clustered at the top right corner vertical wall and parallel to the heated and cooled walls. The heat is transferred to the entire space of the cavity mainly by conduction. The upper part of the enclosure, right below the external flow, has horizontal stratified isotherms which means that the forced flow is stronger and the flow inside the cavity is dominated by heat conduction mechanism. As Ri increases to 0.1 the contour lines deviate from the middle heated wall to the upper part of the cavity due to the dominance of the natural convection effects. In this case the buoyancy forces become stronger and push the heated air further upstream. From Fig. 6 we observe
4802
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
Fig. 7. Flow structures for Re = 100 (top and middle) and Re = 1000 (bottom), (a) Ri = 0.1, (b) Ri = 1, (c) Ri = 10, (d) Ri = 100, (e) Ri = 0.01 and (f) Ri = 0.1.
that the flow exhibits a three-dimensional structure. As the Ri increases further (Ri P 1) the flow becomes unsteady. The temperature profile indicates the Kelvin–Helmholtz instabilities in the shear layer, at Ri = 1, due to the interaction between the forced flow and the buoyancy-induced flow by the natural convection. This instability is located right below the external flow close to the heated wall similar to the one noticed by Yao et al. [18] at high Re number without heat transfer. The flow becomes extremely unstable with unsteady solution. This is why the results are not presented here. Note that the isothermal wall remain parallel in the neighborhood of the cooled wall. At Ri = 10 the isothermlines are very close, the thermal boundary layer from the heated wall becomes thinner, and the hot buoyant fluid rises like a buoyant plume into the outside domain and escape from the restrictive vertical cavity walls. The thickness of this thermal layer decreases as Ri increases, as seen in Fig. 5(c) and (d) as in [13,15]. We found four grid nodes inside this thermal layer. The isotherm pattern indicates that up to a critical value of the Ri number, the buoyancy forces (temperature gradients) exceed the dissipative effects of viscous forces causing wavy pattern of the isotherms at the outflow channel. Its effect is to enhance heat transfer rate of the external flow leaving the channel and generate more unstable wakes. A thermal plume can also be seen to rise into the floor of the cavity, see Fig. 6.
4.2. Development of velocity components To see the impact of the mixed convection on the flow pattern and to give the idea of recirculating and ascending flow structure, the streamlines, and v, w-velocity profiles at specific planes are displayed in Figs. 7–11. In Fig. 7 we present the streamlines map for Re = 100. For Ri = 0.1 in Fig. 7(a), the cold forced flow passes through the cavity without penetrating and creates a recirculating structure that fills the entire cavity. As the Ri was increased, in addition to the primary eddy two other minor eddies can be seen at the left and right bottom corner, see Fig. 7(b), since the heated and cooled airs move upward and downward, respectively, due to the buoyancy forces. For Ri = 10 and 100, a recirculating cell is present at the upper left of the cavity. The temperature gradients becomes stronger and part of the heated flow leaves close to the heated wall and the horizontal adiabatic wall of the channel behind the enclosure. Part of the forced flow also penetrates further downstream in the enclosure. These explain the higher heated transfer already observed in the temperature profiles. Fig. 8 shows the transverse velocity component (v/U0) as a function of the streamwise distance (x/L) at different constant spanwise and vertical plans. Near the side wall (x/L = 0), and for low Ri numbers (see Fig. 8(a) and (b)), the v-velocity component profile presents positive values meaning that the fluid is moving up and
4803
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808 0.15
0.15
z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
0.1
0.1 0.05
v/U
0.05
v/U
z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
0
0
-0.05
-0.05
-0.1
-0.1
-0.15
-0.15 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
x/L 0.3
0.6
0.8
1
x/L z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
1.5
0.2 1
v/U
v/U
0.1
0
0.5
0
-0.1
-0.5
-0.2
-1 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
x/L
0.6
0.8
1
x/L
Fig. 8. v-velocity component for Re = 100 at x- and y-planes, (a) Ri = 0.1, (b) Ri = 1, (c) Ri = 10 and (d) Ri = 100.
0.06
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.04 0.02
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.04 0.02 0
-0.02
v/U
v/U
0
-0.02
-0.04
-0.04
-0.06
-0.06
-0.08
-0.08
-0.1
-0.1 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
z/L 1
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.3
0.6
0.8
1
z/L x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.8 0.6
0.2
v/U
v/U
0.4 0.1
0.2 0
0
-0.2 -0.1 -0.4 -0.2
-0.6 0
0.2
0.4
0.6
z/L
0.8
1
0
0.2
0.4
0.6
z/L
Fig. 9. v-velocity component for Re = 100 at y- and z-planes, (a) Ri = 0.1, (b) Ri = 1, (c) Ri = 10 and (d) Ri = 100.
0.8
1
4804
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
0.06
0.06
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.04
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.04
0.02
w/U
w/U
0.02
0
0
-0.02
-0.02
-0.04
-0.04 0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
z/L 0.2
0.3
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.15
0.6
0.8
1
z/L
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.2
0.1
0.05
w/U
w/U
0.1
0
0
-0.05 -0.1 -0.1 -0.15
-0.2 0
0.2
0.4
0.6
0.8
1
z/L
0
0.2
0.4
0.6
0.8
1
z/L
Fig. 10. w-velocity component for Re = 100 at y- and z-planes, (a) Ri = 0.1, (b) Ri = 1, (c) Ri = 10 and (d) Ri = 100.
negative values near the side wall (x/L = 1) meaning that the flow is moving down up. This leads to the primary recirculation, similar to the natural convection profiles, already observed in the recirculating nature in the flow structure pictures. Near the bottom wall (y/L = 0.25) the convective flow remains almost constant for low Ri values, see Fig. 8(a) and (b). Fig. 8(b) indicates a slight increase in the v-component, for Ri = 0.1, near the side wall (x/L = 1) meaning a quasi-two-dimensional secondary eddy. At Ri = 1 and 10, see Fig. 8(c) and (d), the buoyancy effects considerably modify the v-velocity. The variation is much steeper near the heated wall (x/L = 1). The heated air rises fast with increasing Ri. For Ri = 1, the change of the v-velocity is small at the bottom and middle part near the side wall (x/L = 0), as opposed to the top left side of the cavity, where the v-component experiences a change in the slop. At the bottom left side of the enclosure the v-component presents negative values meaning that the flow is sinking along the cooled wall. At Ri = 10, the vertical velocity is almost constant near the side wall (x/L = 0). The distribution gives the idea of recirculating and ascending flow structure. Similar behaviors were observed in the two-dimensional results of Manza et al. [7] at Ri = 0.1 and Ri = 100 for Re = 100. Fig. 9 shows the v-velocity component as a function of the spanwise distance (z/L) at two different planes, the constant plan x/L = 0.5 at the middle of the cavity and the plan x/L = 0.75 close to the step. The value y/L = 0.75 corresponds to a transverse plan near the top enclosure, while y/L = 0.5 corresponds to a transverse plan at the middle cavity. At the middle constant (x/L = 0.5) plan and for both Ri = 0.01 and 0.1, the v-velocity distribution shows that the flow remains nearly constant. It descends toward the bottom as Ri increases to 1 and 10. Close to the step, the flow ascends toward the top in the central direction for Ri = 0.1 and descends
downward for high Ri values. Near the top enclosure and the heated wall the flow is displaced towards the bottom in the central direction for Ri = 0.01 and 0.1. At Ri = 1 and 10, the vertical velocity component behavior is similar to Fig. 9(a) and (b) discussed above, the flow is ascending in the vicinity of the side walls meaning that the hot air is rising, see also Fig. 6. Fig. 10 shows the w-velocity as a function of the spanwise distance (z/L) at the same constant plans as in Fig. 9. At the middle plan, the flow is directed toward the center of the enclosure. Near the step, the same figure shows that the flow is moving toward the side walls. Even at low Ri, where the flow is quasi-two-dimensional, the movement is always three-dimensional with a non negligible w-component (see Fig. 10(a) and (b)). Flow structures for Ri = 0.01 and 0.1 at Re = 1000 are depicted in Fig. 7(c) and (d). For Ri = 0.01 the streamlines have the same behavior as for Re = 100. The amount of recirculation flow increases due to stronger forced flow effects. The vertical velocity profile in Fig. 11(c) highlights the recirculating cell. For Ri = 0.1 one can observe the development of two recirculating structures. One cell is located just below the airstream and the other cell is in the downstream of the enclosure due to the raised heated air and the sinking cooled air. Both eddies are accompanied with a stratified air at the floor of the cavity just below streamlines of external flow. The velocity component profiles in Fig. 11(d) verify this scenario. Fig. 11 shows also the v- and w-velocity components as a function of the spanwise distance (z/L) and streamwise distance (x/L) at two at two different planes, at the middle and close to the floor cavity. It can be seen that the vertical velocity v, see Fig. 11(a) and (b), increases in magnitude because the forced flow moves faster with increasing Re. The velocities inside the enclosure experiences irregular ascending and descending at the vertical and
4805
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
0.15
0.15
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.1
0.1 0.05
0
v/U
v/U
0.05
-0.05
0 -0.05
-0.1
-0.1
-0.15 -0.2
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0
0.2
0.4
0.6
0.8
-0.15
1
0
0.2
0.4
z/L 0.15
0.2
z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
0.1
1
0.1
0 -0.05
v/U
v/U
0.8
z = 0.5, y = 0.25 z = 0.5, y = 0.5 z = 0.5, y = 0.75
0.15
0.05
-0.1 -0.15
0.05 0 -0.05
-0.2
-0.1
-0.25 -0.3
0.6
z/L
0
0.2
0.4
0.6
0.8
-0.15
1
0
0.2
0.4
x/L 0.06
0.08
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.04
0.6
0.8
1
x/L
x = 0.5, y = 0.5 x = 0.5, y = 0.75 x = 0.75, y = 0.5 x = 0.75, y = 0.75
0.06 0.04
w/U
w/U
0.02
0
0.02 0 -0.02
-0.02
-0.04 -0.04
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
z/L
0.6
0.8
1
z/L
Fig. 11. Velocity components for Re = 1000, left: Ri = 0.01 and right: Ri = 0.1.
horizontal directions and the fluid is no longer steady, see also Fig. 11(f). The w-velocity component, Fig. 11(e), displays the same behavior as in the previous case with Re = 100. The effect of varying the ratio of the height of the vertical heated wall D is shown in Figs. 12(a) and (c) and 13(a) and (c) for Re = 100 and Ri = 0.1. Figs. 12(a) and (c) and 13(a) and (c) show the isotherms and streamlines, respectively, for three cavities of aspect ratio AR = 0.5, 1.5 and 2. For AR = 0.5, the main vortex fills the cavity close to the opening vertical wall. For a deeper cavity with AR = 1.5, the isotherms remains parallel to the heated wall due to the high gradients in the vertical direction. In addition, the intensity of the cavity flow decreases at the bottom wall and the fluid flows along the lower plate. For AR = 2, a pair of counter-rotating vortices appears in the cavity. This can be attributed to the buoyancy forces which opposes the recirculating flow within
the enclosure. The fluid in this case penetrates into the cavity and reverses its direction in the bottom region. The isotherms in the lower part are almost parallel to the vertical heated wall. The isotherms and streamlines for high Richardson number of Ri = 1 and aspect ratio AR = 2 are presented in Figs. 12(d) and 13(d). As the Ri increases, the second vortex which is caused by the strong buoyancy forces becomes much stronger and occupies more than the half of the open cavity. This recirculating flow is responsible for the penetration and inversion of the fluid at the central half portion of the cavity. The time behavior of flow variables at the trailing edge corner for Ri = 1 is shown in Fig. 14. The flow enters the transient regime and becomes unstable. For both high Re and Ri numbers the buoyancy forces exceed the dissipative effects of viscous forces causing the fluid to rise quickly and forming an unsteady convection and fluctu-
4806
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
number at the hot wall is always seen to increase with Ri for each value of Re. Similar behavior were observed in [10] by the author in their bi-dimensional numerical tests. At low values of Ri number and for both Re = 100 and Re = 1000, the Nusselt number increases slightly since the buoyancy effects are weak and the mixed convection transport mechanism is largely dominated by convective effects. Significant increase is found for high values of Ri, up to certain Ri value for each Re number, since the natural convection was the dominant mode of the mixed convection transport. The continuous increase of the heat transfer with increasing Reynolds for each Ri value is also observed. 5. Conclusions
Fig. 12. Temperature distribution for Re = 100 and Ri = 0.1, (a) D = 0.5, (b) D = 1.5, (c) D = 2, and (d) D = 2 with Ri = 1.
ation in velocity and temperature. This behavior will need more detailed study and would constitute a subject of a future work. 4.3. Heat transfer The averaged Nusselt number calculated at the hot and cooled wall is presented against Ri for Re = 100 and 1000 in Fig. 15. At cooled wall, the Nusslet number decreases slightly with Ri. It becomes very small for Re = 100 and Ri = 100 since the isotherm lines are clustered parallel to the hot wall, see also Fig. 4. The Nusselt
A detailed numerical study has been carried out to investigate mixed convection for incompressible laminar flow past open cubical cavity. The results show that the flow exhibits a three-dimensional structure and becomes steady for Re = 100 with Ri 6 100 and Re = 1000 with Ri 6 0.1. The forced flow dominates the flow transport mechanism and large recirculating zone form inside the enclosure which results in heat transfer by conduction. In this case the Nusselt number increases slightly. For both high Reynolds and Richardson numbers the natural convection comes into play and push the recirculating zone further upstream. The flow becomes unsteady with Kelvin–Helmholtz instabilities at the shear layer and the Nusselt number increases significantly. For Re = 1000 and Gr = 107 the isotherm lines remain parallel close to the heated wall and form and thermal boundary layer. The work also gave quantitative results for the velocity components resulting from various locations to give the idea of recirculating and ascending flow structure. This flow has revealed to be a very convenient case for studying flow transients under the combined effects of forced and natural convection with a progressive evolution from a quasi-two-dimensional flow to a three-dimensional steady and
Fig. 13. Flow structures for Re = 100 and Ri = 0.1, (a) D = 0.5, (b) D = 1.5, (c) D = 2 and (d) D = 2 with Ri = 1.
4807
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
0.6
0.6 0.5
0.5
0.4 0.4 0.3 0.2
V
U
0.3 0.2
0.1 0
0.1
-0.1 0 -0.2 -0.1 -0.2
-0.3 20
40
60
80
100
120
140
160
-0.4
20
40
60
Time
80
100
120
140
160
Time
0.3
0.5
0.2
0.4
0.1
T
W
0.3 0
0.2 -0.1
0.1
-0.2
-0.3
20
40
60
80
100
120
140
160
0
20
40
Time
60
80
100
120
140
160
Time
Fig. 14. Time history of the instantaneous velocity components and temperature for Re = 1000 and Ri = 1.
16
Acknowledgments
14
This work was financially supported by project DPI2006-02477 from Ministerio de Educación y Ciencia (Spain) and Fondos FEDER.
12 10
References
Nu
8 6 4 2 0 -2 -4 10
-4
10
-3
10
-2
10
-1
10
0
10
1
10
2
10
3
Ri Fig. 15. Average Nusselt number vs. Ri numbers. Cold wall, ‘D’: Re = 100 and ‘r’: Re = 1000. Hot wall, ‘’: Re = 100 and ‘}’: Re = 1000.
unsteady flow in a relatively tidy range of Re and Ri numbers. Future work are being undertaken to study the unsteady behavior for both high Re and Ri.
[1] E. Bilgen, H. Oztop, Natural convection heat transfer in partially open inclined square cavities, Int. J. Heat Mass Transfer 48 (2005) 1470–1479. [2] E. Chakroun, M.M. Elsayed, S.F. Al-Fahed, Experimental measurements of heat transfer coefficient in partially/fully opened titled partially open cavities, J. Solar Eng. 119 (1997) 298–303. [3] Y.L. Chan, C.L. Tien, A numerical study of two dimensional natural convection in square cavities, Numer. Heat Transfer 8 (1985) 65–80. [4] P. Gaskell, A. Lau, Curvature-compensated convective transport: SMART, a new boundedness-preserving transport algorithm, Int. J. Numer. Meth. Fluids 8 (1988) 617–641. [5] A. Javam, S.W. Armfield, Stability and transition of stratified natural convection flow in open cavities, J. Fluid Mech. 445 (2001) 285–303. [6] K. Khanafer, K. Vafai, M. Lightstone, Mixed convection heat transfer in twodimensional open-ended enclosures, Int. J. Heat Mass Transfer 45 (2002) 5171–5190. [7] O. Manza, S. Nardini, K. Khanafer, K. Vafai, Effect of heated wall position on mixed convection with an open cavity, Numer. Heat Transfer, Part A 43 (2003) 259–282. [8] O. Manza, S. Nardini, K. Vafai, Experimental investigation of mixed convection in a channel with an open cavity, Exp. Heat Transfer 19 (2006) 53–68. [9] J. Pallares, I. Cuesta, F.X. Grau, F. Giralt, Natural convection in a cubical cavity heated from bellow at low Rayleigh numbers, Int. J. Heat Mass Transfer 39 (1996) 3233–3247.
4808
Y. Stiriba et al. / International Journal of Heat and Mass Transfer 53 (2010) 4797–4808
[10] E. Papanicolaou, Y. Jaluria, Mixed convection from an isolated heat source in a rectangular enclosure, Numer. Heat Transfer, Part A: Appl. 18 (4) (1990) 427– 461. [11] E. Papanicolaou, Y. Jaluria, Transition to a periodic regime in mixed convection in a square cavity, J. Fluid Mech. 239 (1992) 489–509. [12] E. Papanicolaou, Y. Jaluria, Mixed convection from an localized heat source in a cavity with conducting walls, Numer. Heat Transfer, Part A: Applications 23 (4) (1990) 463–484. [13] K. Vafai, J. Ettefagh, The effects of sharp corners on buoyancy-driven flow with particular emphasis on outer boundaries, Int. J. Heat Mass Transfer 33 (1990) 2318–2328. [14] K. Vafai, J. Ettefagh, Thermal and fluid flow instabilities in buoyancy-driven flows in open-ended cavities, Int. J. Heat Mass Transfer 33 (1990) 2329–2344.
[15] K. Vafai, J. Ettefagh, Thermal and fluid flow instabilities in buoyancy-driven flows in open-ended cavities, ASME J. Heat Transfer 116 (1994) 604–613. [16] K. Vafai, P.C. Huang, Analysis of heat transfer regulation and modification employing intermittency emplaced porous cavities, Int. J. Heat Mass Transfer 33 (1990) 2329–2344. [17] H.A. Van der Vost, BiCGSTAB: a fast and smooth converging variant of Bi-CG for the solution of non-symmetric linear system, SIAM J. Sci. Statist. Comput. 13 (1992) 631–644. [18] H. Yao, R.K. Cooper, S. Raghunathan, Numerical simulation of incompressible laminar flow over three-dimensional rectangular cavities, J. Fluid Eng. 126 (2004) 919. [19] P.S.B. Zdanski, M.A. Ortega Jr., N.G.C.R. Fico, Heat transfer in the flow over shallow cavities, J. Heat Transfer ASME 127 (2005) 699–712.