1. Introduction
From the pioneering experimental study of Elenbaas [
1], configurations formed by vertical heated plates with fluid flows induced by buoyancy effects have been the topic of several studies (Incropera and De Witt [
2], Bejan [
3]). The natural convection is a mode of heat transfer that presents undoubted advantages, under given circumstances. The vertical channel configuration with isothermal heating (
Figure 1a) can be regarded as the most typical analyzed configuration. A relevant review of benchmark solutions for natural convection flows in vertical channels was presented recently by Desrayaud et al. [
4].
The
optimization problem for the described issue is represented by the electronic equipment cooling. In this case, one of the most important objectives consists of dissipating the heat generated in the devices to avoid overheating. The determination of the thermally optimum spacing
bopt (between the plates or walls forming the vertical channel) that maximizes the heat transfer rate per unit area can achieve the
thermal optimization of the system. The first criterion for obtaining
bopt was reported by Bodoia and Osterle [
5]. Analytically, a given
optimization function can be obtained from the appropriate correlations for calculating the heat transfer rate. Bar-Cohen and Rohsenow [
6] and Zamora and Hernández [
7] employed different correlations for the average Nusselt number for achieving this aim, with laminar flow assumption. The optimum was reached for a modified Rayleigh number Ra* (based on the inter-plate spacing) placed near the crossing point between the fully-developed and the boundary layer asymptotes; in this manner, its value depends on the fitting constants used in the average Nusselt number correlations. In reality, the idea of considering an isothermal vertical channel as a configuration for obtaining the optimal spacing between plates makes physical sense if this geometry is regarded as a representative sample of the channels formed with equally-spaced isothermal plates located into an available horizontal gap
L. When the number of plates increases, the total area for heat transmission increases, but in turn the temperature gradient at the walls decreases, as well as the average heat transfer coefficient. Since the described effects become the opposite when the number of plates decreases, then a thermally optimum wall-to-wall spacing can be encountered. Note that if the
theoretical bopt is
bopt <
L, a sole channel must be formed. The global heat transfer could not be maximum in this case, but anyhow its value will be higher than that reached generating more than one channel in the available space.
In some cases (in passive ventilation systems, for instance), the
dynamic optimization is clearly interesting. Here, obtaining the maximum induced mass-flow rate can be the main objective for ventilation purposes (Zamora and Kaiser [
8], Zavala-Guillén et al. [
9]). Previously, the authors carried out calculations for obtaining both thermal and dynamic optimization in geometries concerned to thermal passive systems of buildings (see Zamora and Kaiser [
8]).
Other aspects of the problem can be treated; for instance, finned surfaces in natural convection provide relevant improvements in the thermal performance. Therefore, the distribution of heat sources can be analyzed for achieving the optimal arrangement (da Silva et al. [
10], Zhang and Liu [
11]).
More realistic configurations should be regarded for obtaining practical results. As expected, several works dealing with the geometric optimization of cavity configurations for obtaining the maximum heat transfer by buoyancy effects can be found in the literature (Aounallah et al. [
12], for cavities with curved walls; Biserni et al. [
13], for H-shaped cavities; Lorenzini et al. [
14], for T-shaped cavities; Lorenzini and Rocha [
15], for cavities with T-Y shape, among others). A feasible analysis procedure consists of using both the
scale analysis method and the intersection of the asymptotes obtained for different regimes of the flow (da Silva and Gosselin [
16], for instance). This way is similar to that followed by Bar-Cohen and Rohsenow [
6] and Zamora and Hernández [
7]. Another method is based on the
Constructal Law of Bejan (Bejan [
17], Bejan and Lorente [
18]). These authors pointed out that, essentially, the constructal theory consists of global objectives and global constraints, whereas the geometry of the flow is unknown, i.e., the geometry is not previously assumed, but it is deduced. Alternative configurations and shapes are considered in the literature (da Silva and Gosselin [
16], among others). Two configurations for the buoyancy-induced airflows in vertical channels have been studied in a previous work (Zamora [
19]); one a smooth, curved channel, and the other a channel with several sharp changes of direction. For any given circumstances, the performance of the alternative outlined configurations can be considered better than that corresponding to a straight channel. This is a representative case, in which the geometries are mainly defined, and generalized enough correlations for the Nusselt number are not available; therefore, massive computational data (or experimental, if applicable) should be analyzed in order to find possible optimal situations.
For simulating the buoyancy effects, the well-known Boussinesq approach assumes constant properties of fluid, except for density variations exclusively due to temperature variations in the vertical component of the momentum equation. When the influence of variation of the
air thermophysical properties is taken into account, a contrasted result found in the literature is that heat transfer coefficients and the induced mass-flow rate are considerably lower than those obtained under non-Boussinesq conditions. This particular effect can be attributed to the increase of the air viscosity (
viscous drag) and the decrease of the air density, which produces an additional thermally-induced pressure drop (
thermal drag), as the temperature difference increases (Zhong et al. [
20], Guo and Wu [
21]). In other words, both the viscous drag and the thermal drag increase faster than the buoyancy force when the temperature increases. The flow patterns could be affected strongly by the fluid variable properties for intense heating conditions; in this way, it can be expected that the thermally optimum inter-plate spacing (if it exists) is also affected by the described phenomenon.
In present work, a realistic situation that consists of considering a limited space, i.e., a vented cavity, into which a vertical channel is placed, is outlined (see
Figure 1b). The spacing between the plates forming the channel can vary, but the dimensions of the cavity remain constant. The physically periodic conditions of the array of plates described above disappear, and thus the geometrical surroundings of the isothermal channel have a relevant influence on the behavior of the convective airflow. In this case, the area dedicated to convective heat transfer does not vary, but in turn, the thickness of the thermal boundary layers adjacent to heated walls varies strongly when the inter-plate spacing
b changes. Therefore, a thermally optimum wall-to-wall spacing could exist; its determination constitutes the aim motivation of this work, as well as the simulation of the evolution of air properties. The adopted approach is numerical; results are presented as a function of the inter-plate spacing for a wide range of Rayleigh numbers. In view of the fact that the main applications of the outlined configuration might be subject to intense heating conditions, the influence of considering variable thermophysical properties of air is also included in the study. It is foreseeable that some effects could be explained by the flow configuration changes when the relevant parameters vary.
2. Materials and Methods
Figure 1 shows the
configurations outlined in this paper. The first one is a typical vertical channel (
Figure 1a); the height is
H (
y-direction) and the spacing between walls is
b (
x-direction); the relationship
b/
H is commonly named the
aspect ratio.
Figure 1b shows the configuration proposed for the numerical simulations; the dimensions of the cavity are
HT in
y-direction and
L in
x-direction, whereas the vents have a width equal to
bV; a standard vertical channel (
b ×
H) is placed centrally into the cavity; the thickness of the plates is considered negligible. A square cavity,
L/
HT = 1, is considered; the ratio between channel and cavity heights is
H/HT = 1/2, and the width of vents is given by
bV/
HT = 0.075. The flow enters mainly through the lower vents and goes out mainly through the upper one.
The plates of the channel are considered isothermal at Tw, whereas the reference (ambient) conditions for air are given by p∞ = 105 N/m2 and T∞ = 293 K. In this manner, the relevant parameters to be considered are based on the height of plates H and on the characteristic temperature difference ΔT = Tw − T∞. Hence, the Rayleigh number is defined as RaH = GrH Pr∞, being GrH the Grashof number based on H, GrH = gβ (Tw − T∞) H3/ν∞2, with g the gravity acceleration, β = 1/T∞ (perfect gas assumption), ν∞ the kinematic viscosity of the fluid, and Pr the Prandtl number of the fluid (Pr∞ = 0.7 in this work). The range of the Rayleigh number considered is 103–1011 with turbulent simulation, which includes the whole range of typical performance conditions.
The
average Nusselt number at heated, isothermal walls is calculated as follows
The Nusselt number calculated at the plate side facing the surrounding cavity wall is denoted with subscript out, and that calculated at the side facing the other plate is denoted with in.
The buoyancy-induced incoming (or outgoing) dimensionless
mass-flow rate M is
m being the two-dimensional mass-flow rate.
Lastly, for evaluating the fluid variable properties effects, the
heating parameter can be defined as follows
For low enough values of Λ, the influence of variable properties is negligible, and therefore the Boussinesq approach can be employed. In this work, this limit is achieved for the lowest value considered of the heating parameter, Λ = 0.01. For analyzing the variable properties effects, the range of values of Λ considered for numerical simulations is 0.01–10, which includes the whole range of typical situations.
Regarding the
mathematical approach, the two-dimensional form of the elliptic time-averaged Navier–Stokes equations (continuity, momentum, and energy) for the turbulent airflow (see Zamora and Kaiser [
8] for more details), are solved numerically by using the general-purpose Phoenics code (Phoenics Encyclopaedia [
22]), which is based on a finite-volume procedure. The
driving (
buoyancy) force is calculated directly from density variations (instead of from temperature variations such as is advocated by the Boussinesq approach), and it is implemented in the vertical component of the momentum equation. The selected turbulence model for solving the closure problem is the version of the two-transport equations
k–
ω model proposed by Wilcox [
23], which includes a low-Reynolds extension for near-wall turbulence. If the mesh includes some nodes into the boundary layer, the entire laminar–transitional–turbulent behavior can be appropriately simulated. The obtained results have been compared with those computed through the SST
k–
ω turbulence model, but they are not presented because of differences tend to be negligible provided that the grid was refined enough, mainly near the walls. Now then, the necessary computations times are, in general, higher with the SST
k–
ω model.
Regarding the boundary conditions, the plates of the channel are considered with a uniform hot temperature T = Tw, whereas a cold temperature T = T∞ is imposed on the fluid outside the cavity. The walls forming the cavity are treated as adiabatic. The non-slip condition for velocity and turbulent kinetic energy is imposed at all the walls. At the inlet sections (lower vents), the mass-flow rate is assumed to be dependent of the square root of the difference between the ambient pressure p∞ and the pressure p computed at each inlet cell (this can be regarded as a result of the Bernoulli equation application at the surroundings of the cavity), whereas the temperature is fixed at T∞. At the outlet sections (upper vents), the pressure is fixed at p∞, and the streamwise variations of velocity components, temperature, and turbulent variables are neglected.
From the experimental data available in the literature (Çengel and Cimbala [
24], for instance), own correlations are developed for expressing the evolution of the thermophysical properties of air as a function of temperature. Thermal conductivity, specific heat at constant pressure and kinematic viscosity are increasing with absolute temperature
T,
valid for 273–2273 K, with determination coefficients R
2 equal to 0.9937, 0.9993, and 0.9990, respectively. In addition, the perfect gas assumption
ρ =
p/
RT, being
R the air constant, equal to 287 J/kg.K indicates that since variations of pressure are low, the density of air is mainly decreasing with
T. These correlations are implemented into the Phoenics code for simulating the thermophysical air behavior.
Regarding the numerical details, the time-dependent procedure is employed for obtaining the numerical results; a fully-implicit scheme for the time discretization is used. A scale analysis in the y-momentum equation gives a typical vertical velocity in the order V ≈ [gβ(ΔT)H]1/2 for natural convection flows, by equating convective and buoyancy terms, v(∂v/∂y) ≈ gβ(ΔT). The characteristic time can be estimated as t0 ≈ H/V. For obtaining the appropriate accuracy, the employed values of time step Δt are into the range 0.01–0.25. A number of 40 iterations allows us to achieve spatial convergence in each time step. Finally, most of the computations lead to obtaining essentially steady-state solutions.
The equations are discretized by a staggered-grid scheme. Structured, cartesian meshes are employed (120 × 120 cells in most cases), with different power-law distributions to produce fine meshing near the walls, as aforementioned. A typical grid is shown in
Figure 1c. The accuracy of the numerical results is tested by a grid dependence study; both the dimensionless sub-layer scaled distance
y+ and the total number of cells are taken into account for analyzing the influence of meshing on the results (see Zamora and Kaiser [
8]). Parameter
y+ is revealed as the most influential so that for low enough values of
y+, the obtained results are affected mainly by the given value of
y+. Therefore, the results presented are computed using grids with
y+ low enough (mainly in the range 0.1–0.5), achieving in this way the grid independence.
The numerical results are obtained through a second-order “muscl” (non-linear) differencing scheme (Phoenics Encyclopaedia [
22], Van Leer [
25]) for the convective terms of the momentum, energy, and turbulence transport equations. Focusing on numerical convergence, the relative change of any dependent variable in each iteration is less than 10
−5. The relative residuals for mass, momentum, energy, and turbulent variables for the full flow field are less than 10
−4.
4. Discussion of Flow Patterns
Flow patterns are shown in
Figure 9,
Figure 10 and
Figure 11.
Figure 9 illustrates the thermal flow patterns that are found for two different Rayleigh numbers and two different aspect ratios. Here, it can be considered that results are obtained as the Boussinesq approach proposes since that Λ = 0.01, and therefore the air properties are considered almost constant. Note that for low enough values of Ra
H, the temperature field tends to be more homogeneous throughout the cavity than for higher values. The buoyant plumes generated at the upper part of heated plates are clearly more recognizable when
b/
H is low enough, and for
b/
H high enough, plumes are rapidly routed towards the upper exit vents. Characteristic values of dimensionless temperature difference
θ = (
T −
T∞)/(
Tw −
T∞) are depicted in this figure.
The effects of variable thermophysical properties are graphically depicted in
Figure 10 and
Figure 11. The first shows the filled contours of temperature, and the second, streamlines and selected velocity profiles, for a characteristic value of Rayleigh number Ra
H =10
9, aspect ratio
b/
H = 0.5, and different values of heating parameter Λ. In
Figure 10, it is evident the progressive thickening of the thermal boundary layers adjacent to heated walls as Λ increases. The isotherms tend to occupy a larger portion of the cavity for higher values of Λ, and the Nusselt numbers at walls tend to decrease as a consequence of the generated thermal drag for intense heating conditions. This behavior is accompanied by a relevant change in the velocity field of airflow, which can be observed in
Figure 11. When Λ increases, the increase of viscosity produces a reduction of the gradients of velocity into the boundary layers adjacent to walls; in addition, the generated viscous drag, as well as the decrease of buoyancy forces, modify the global field of velocity, and therefore the velocity profiles through the immersed channel. Both the average velocity and the density decrease significantly for intense heating conditions, and, finally, the induced mass-flow rate decreases. In this figure, note that since reached values of velocity are very low in the upper part of the cavity, then some streamlines seem to end on the top wall. In
Figure 11, the maximum values of the dimensionless vertical velocity ψ =
VH/
νGr
H, computed in each case, are indicated.
In general, the structure of the airflow consists of two clearly defined incoming jets that rise through the core of cavity, impelled by the buoyancy forces produced by the heated immersed channel (see
Figure 11). Two counterrotating rolls are detected in most of the cases, placed at the spaces comprised between the outer sides of channel plates and the walls of the cavity. Under given circumstances, some additional recirculation areas are found at the central, upper, and bottom parts of the cavity. An asymmetric alternative configuration of the flow motion is encountered for some intermediate values of both Ra
H and (not too high)
b/
H. This flow pattern consists of a predominant jet from the right that crosses the cavity and goes out through the upper exit vent placed on the left (or vice versa). There are well-founded suspicions that this phenomenon could be a numerical effect; in fact, reducing the time step of numerical marching, the explained asymmetry tends to be reduced as well; however, it can be noted that these oscillating flows could reproduce a certain physical reality, which should be analyzed in detail in each case.