Numerical investigation of free convection through a horizontal open-ended axisymmetric cavity

Objectives : The purposes of this work are to investigate the free convective heat transfer in an axis-symmetric open-ended cavity heated from below and to propose useful correlations of Nusselt number. Methods : The governing equations that model the ﬂuid ﬂow and the temperature ﬁeld are solved using a control volume-based ﬁnite diﬀerences method. Under steady state condition, the natural convective ﬂow is considered to be laminar, incompressible and axisymmetric. The Boussinesq approximation with constant thermophysical properties is adopted. Numerical experimentations are performed to deduce the optimum sizes of the calculation domain and the mesh grid. Findings: the obtained results indicate that when Rayleigh number (Ra) and aspect ratio (A) are low the heat transfer is weak and mainly conductive. The increase of Ra and A enhances the convective heat transfer mode thereby the heat transfer is ameliorated. Unlike the Rayleigh Bénard convection, the transition from conduction to convection produces at critical value of Rayleigh number (Ra c ) that is variable dependent on A. Novelty: To the best of authors knowledge, the formula of (Ra c ) elaborated in this work for the studied cavity is the ﬁrst attempt. As well, correlation of Nusselt numbers (Nu) for the cold upper plate in terms of Ra and A is performed. Comparisons between Nu at the lower plate given in previous work and Nusselt number at the upper plate is conducted.


Introduction
The natural convection heat transfer plays an important role in many applications and industrial processes such as electronics equipment, heat exchangers, and Chemical Vapor Deposition (CVD) reactors. In fact, the heat dissipation of electronics equipment is a major industrial interest seeing the trend toward miniaturization of equipment. Certainly, the use of forced convection by means of fans is very efficient. But this approach requires more space with an increase in energy consumption and noise generation. The fundamental knowledge of the free convection flow and the characteristics of heat transfer are important for innovation in practical applications. However, several challenges are imposed in its study due to the high thermal and hydro-dynamical coupling as well as the strong sensitivity to external disturbances. These difficulties grow when flows are encountered in complex geometries such as open cavities. To date, the majority of the works conducted concern the rectangular cavities (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12) . Experimental flow visualization and air temperature measurements were employed by Manca and Nardini (5) to describe the natural convection inside an open-ended cavity with a heated upper plate and an unheated lower plate. The same configuration was investigated numerically by Andreozi and Manca (6) using Boussinesq approximation under steady-state. Their numerical simulations illustrated that the fluid penetrated inside the cavity up to the midplane with a C-loop pattern. Boetcher and Sparrow (7) accomplished a numerical study of buoyancy induced flow in a horizontal open-ended cavity by means of ANSYS CFX 11.0 software. They discussed thoroughly the relevant numerical issues using enlarged calculation domain, appropriate boundary conditions on the surfaces of the extended domain and meshing to achieve high accuracy. Kitamura and Asakawa (10) have conducted an experimental study in a water medium heated from below. They have brought to light thermal instabilities that could occur. Fu et al. (11) conducted a numerical investigation of natural convection in three-dimensional horizontal parallel plates with a heated bottom surface without using Boussinesq approximation. The results show that both thermal and flow fields become unsteady when Rayleigh number growth closely to 3.15 × 10 6 . In another numerical study, Fu et al. (12) concluded that more the space between the plates is narrow more the thermal field and the fluid flow are stable.
On the contrary, the axis-symmetric cavities have attracted less interest especially the open ones. However, these cavities are widely present in many electronics equipment such as power component with cylindrical shape facing a casing wall. Van Santen et al. (13) conducted a numerical study of the onset of thermal instability as a function of Reynolds number, Prandtl (Pr) and Rayleigh (Ra) between two circular plates differentially heated. The study was conducted on a mixed convection flow. An experimental study was conducted by Hideki et al. (14) for natural convection heat transfer in a narrow gap between two horizontal circular plates in air heated from below. The heat transfer measurement and the flow visualization found to be in a good agreement with the results of the numerical study conducted later by Hideki et al (15) of the same configuration. Horibe et al. (16) presented an experimental study of natural convection in an enclosure containing two horizontal parallel heated circular plates. The authors had shown that the flow and the heat transfer are strongly affected by the space between the plates and by the heat flux. Chehouani et al. (17) performed an experimental study to investigate heat transfer in an open ended cavity between two horizontal parallel circular plates. The lower plate is isothermally heated whereas the upper plate is maintained at ambient temperature. The same configuration studied in (17) was investigated both numerically and experimentally using holographic interferometry by AIT HAJ SAID et al. (18) , qualitative similarity between the numerical and experimental results was observed and the increase of the heat transfer is associated with the increase of the Rayleigh number and aspect ratio. Yanping et al. (19) presented a numerical investigation on the natural convection heat transfer in a circular enclosure with an internal cylinder; the synergy analysis conducted in this work reveals that the heat transfer depends significantly on the form and orientations of the cylinder. More recently Medebber et al (20) conducted a numerical investigation of transient free convection in a vertical cylinder partially annulus, uniform temperature is imposed cross vertical walls, while the top and bottom walls are adiabatic. In this work the authors concluded that the heat transfer is great during the first moment. When the steady state is reached, the heat transfer is reduced.
Most of the works cited above was interested by giving useful correlations of Nusselt number. Table 1 summarizes these correlations with their conditions of elaboration and domain of validation.
In this work, we performed a numerical study of the convective heat transfer in an axisymmetric cavity with vertical gradient of temperature. The considered system is that of two circular plates suspended in an infinite space and maintained at fixed temperatures. The lower plate is heated at uniformly constant temperature and the upper plate is kept at ambient temperature. The objectives of our investigation are to bring up and explain the effect of Rayleigh number and aspect ratio on heat transfer characteristics and propose useful correlations between Nusselt number, Rayleigh number and aspect ratio.

Problem definition and mathematical formulation
The system analyzed is composed of two horizontal parallel circular plates, with the upper plate maintained at T 0 , and the lower plate is heated at T H with T H > T 0 . As shown in Figure 1, H is the distance between the two plates with the same radius R.
https://www.indjst.org/ The natural convective flow between the two plates is considered to be laminar, incompressible and 2D-axisymmetric. Under steady state condition, the governing equations in the non-dimensional form, adopting Boussinesq approximation with constant thermophysical properties, are: Continuity equation: Motion equations: Radial component Axial component Energy equation The governing equations (1) -(4) are made dimensionless by scaling the lengths variables by H, the velocities by α H and the pressure by ρ α 2 H 2 .
Defining the stream function and vorticity (Ψ, ω) as follows: Instead of equations (1) -(4), the equations that have to be solved are: The local Nusselt number is defined as follows h is the local heat transfer coefficient obtained through energy balance on a surface as : Where T s is the surface temperature, T r is the reference temperature and n denotes the outward pointing normal from the surface over which the heat flux is calculated. The mean bulk temperature is usually taken as the reference temperature T r . In the problem considered here the mean bulk temperature approaches the ambient temperature T 0 . As a result, for the hot surface, Eq. (11) can be written as: The same equation cannot be written for the cold surface as it is set to the same temperature as the ambient temperature. As a result, the heat transfer coefficients for the hot and cold surfaces are calculated using the temperature difference ∆T between the hot and cold surfaces. Introducing the non-dimensional temperature, we obtain the following formula to compute the heat transfer coefficient: h is positive for the lower plate and negative for the upper one. This justifies the added symbol of absolute value in Eq. (10). The mean Nusselt numbers for lower circular plate Nu(low) and upper circular plate Nu(up) are given by: This number is variable depending on Rayleigh number, Prandtl number and the aspect ratio.

Numerical resolution method
The governing equations in terms of (ψ, ω, θ ) in the mathematical model are of the convection-diffusion type and were written in the general form by Gosman (21) . This common formulation has many practical advantages since it allows the use of one numerical algorithm to solve all equations. The computational domain is discretized into a number of adjoining non-uniform https://www.indjst.org/ cells. The model equations are solved using a control volume-based finite differences method described in (21) . The general equation is integrated over the surface element surrounding each point of the mesh. The sets of algebraic equations resulting from the numerical scheme were solved iteratively using the Gauss-Seidel method.
It is assumed that the solution has been reached when the following convergence criteria ϕ l+1 −ϕ l ϕ l ≤ 10 −6 for temperature, stream function and vorticity have been met. In the previous expression, l refers to a number of iterations and ϕ represents one of the three dependent variables.
The mathematical model and the numerical resolution described above led to the elaboration of a C++ computer code which was successfully and well validated in previous work (18) .

Calculation domain and boundary conditions 4.1 Calculation domain
In the numerical analysis of the natural convection in open-ended cavity problems, the major difficulty lies in specifying the adequate boundary conditions at the cavity aperture. In this study, we preferred to use extension domain rather than another numerical approach (22) . However, this option brings up other quandaries that are, the location of the extended boundaries and the appropriate boundaries conditions to be taken to reflect the physical reality. The commonly adopted approach is to perform a numerical experimentation to deduce the minimum size above which the domain extension has little influence on the solution.

Boundary conditions
The choice of the conditions at the extended boundaries remains a controversy problem. In this study, we assumed that the cavity is placed in a large room open from the top with solid walls (see Figure 1). This assumption leads to a simple mathematical formulation of the boundary conditions and easy numerical implementation using the formalism (ψ, ω) with the drawback that the computational domain must be sufficiently large.
The assumption of axis-symmetry makes it possible to restrict computational domain to half field which reduces the necessary memory and time of calculation to the half. The required size of the extended domain will be discussed later in this paper. The limits of computational domain, where the boundaries conditions are imposed, are presented in Figure 1.
No-slip and impermeable walls, no radial gradient at the centerline and no axial gradient at the outflows are assumed. The plates are assumed to be isothermal. The right and the bottom boundaries of the computational domain correspond to free surfaces at ambient temperature which are located sufficiently far from the cavity axis so that the propagation of boundaryinduced disturbances is minimized. Finally, the horizontal plates are assumed to have a non-null height which is considered adiabatic. Table 1 summarizes all the boundary conditions used in this investigation.

Effect of calculation domain size
The required dimensionless lengths with respect to R, L 1 in r direction and L 2 in z direction, are unknown. They should be long enough so that the heat transfer and the flow become insensitive to their values. Several numerical runs were performed with various domain length L 1 and height L 2 leading to find the adequate domain size which permits a good compromise between the accuracy of the obtained results and the computational cost. Examples of the performed tests are given in Table 2 for A=0.5 https://www.indjst.org/ and different Rayleigh number. Uniform mesh grids with a spacing of 0.05 were used to perform these tests. Comparisons are presented in terms of the sensitivity of − Nu(low), − Nu(up), ψ min and ψ max The parameters ψ max and ψ min are respectively the minimum and maximum of the stream function in the cavity. It was found that the optimum size (L 1 xL 2 ) is (5x7). Indeed, the passage from (2x3) to (5x7) leads to significant variations, about 8.5% for − Nu(low) and 19% for ψ max . Beyond the domain size (5x7), the calculated parameters vary by less than 0.4%. The same tests were performed for aspect ratios A =1 and 0.25. It was carefully verified that (5x7) is sufficient to model the flow and heat transfer for all parameters considered of Ra and A.

Effect of mesh grid size
To verify the effect of grid size, various numerical calculations were done by increasing the number of nodes. A non-uniform mesh grid is used to save computational time. The calculation domain contains two parts; the cavity interior, which is the space between the circular plates, and the cavity exterior which is defined by the whole extended computational domain.
Careful attentions are given to the gridding structure. In the cavity interior, which is the essential part of interest, the mesh grid is very fine near the circular plates because of expected high temperature gradients. Along the plate surfaces, a parabolic mesh grid structure is adopted. So, the grid is fine near to the centerline and in the apertures of the cavity that way, the effect of the sharp corners is overcome (4) . In the cavity exterior, the grid structure is smooth in the vicinity of the cavity and bit by bit it becomes rough toward the far field limit. Extensive tests were conducted for five mesh grid sizes namely M 1 , M 2 , M 3 , M 4 and M 5 . Table 3 shows the characteristics of each mesh grid.  Table 4 represents the results obtained for the considered parameters of Rayleigh numbers under A=0.5. The number of nodes in Table 5 corresponds to the mesh grid of the cavity interior. In Table 5, we present the higher deviations recorded in Table 4.
From Table 5, we stipulate that satisfactory results can be reached using the grid M 4 . Indeed, from M 4 to M 5 the calculated parameters vary by about 0.6% in the worst case for ψ max at Ra=10 5 whereas we note less than 0.12% and almost no variations for Ra=10 4 and Ra=10 3 respectively (Table 4). https://www.indjst.org/

Temperature field and flow structure
In order to be clear and more understandable in this article, we represent a brief description of the flow patterns and the temperature fields at Ra = 10 3 and 10 5 in Figures 2 and 3 for A= 0.2 and 0.5. For more information the reader may consult the previous work (18) where the topology of the medium and the influence of Rayleigh number and aspect ratio were deeply discussed. In Figures 2 and 3, only the space between the circular plates and their vicinity are shown. The rest of domain has been omitted.
https://www.indjst.org/  For weak Ra and A≤0.3, it is noticed that the medium is stratified (Figure 2) and the distribution of the temperature field is horizontal. This indicates that heat transfer is dominated by conduction. As Ra increases, the convective flow becomes more important. Indeed, for A = 0.2 and Ra = 10 5 , the isotherms become more curved due to the appearance of two rotating air areas.
https://www.indjst.org/ For aspect ratios sufficiently high A=0.35 to 1, When Ra=10 3 , the convection is very weak. A few amounts of cold fluid is driven into the cavity. As Ra increases the convective flow is enhanced. The temperature field is distributed vertically throughout the entire cavity except in the vicinity of the two plates where the distribution is horizontal. Thus, a strong temperature gradient is localized at the edge of the lower plate and in the middle of the upper plate.
Qualitatively the observed medium topology in this work is similar to the observed one in (4) where the studied cavity is rectangular, differentially heated and the flow field is assumed to be 2D. In case of high aspect ratio (A ≥ 0.35) similar structure of flow and temperature fields are also obtained in (5) and (6) where the studied cavity is rectangular with adiabatic lower plate and heated upper plate, the numerical study (6) is also performed with the assumption of 2D flow. Therefore it will be interesting to mention that in the current work the convective filed is considered axisymetric as it was experimentally verified in (17) and (18) . So the whole 3D field is obtained by performing the numerical resolution in only 2D.

Local heat transfer
The radial distribution of the local Nusselt number at the surface of the hot bottom plate for different Rayleigh numbers as a function of aspect ratio are shown in Figures 4, 5, 6 and 7.   As expected, the local Nusselt number and its maximum which is located at the edge of the plate in all cases increases with the Rayleigh number. However, the influence of the aspect ratio deserves more discussion. At this level, we notice the effect of three processes. The first is the dominance of the conductive regime in the case of low Rayleigh numbers or for the low aspect ratios in the case of moderate Rayleigh numbers. In this case, the Nusselt number is constant throughout the surface except at the edge of the plate where it takes a maximum value and it increases by decreasing the distance between the two plates. Indeed, the streamlines presented in Figures 2 and 3 show that the incoming cold fluid that penetrates the cavity meets at first the edge of the lower hot plate. The second process is the heat transfer by convection due to incoming cold fluid from outside the cavity. This process is produced for large aspect ratios and also in the immediate opening of the cavity for any aspect ratio. The third and final process is convection due to secondary flow created in the cases of low aspect ratios under moderate to high Rayleigh condition. At this level, the variation of the Nusselt number is periodic or pseudo-periodic in the radial direction with low values of Nu corresponding to the part of a vortex where the fluid is rising and vice versa for higher values of Nu. The amplitude of these oscillations increases while Rayleigh number increases.
https://www.indjst.org/ Finally, we notice that the maximum values of Nusselt number calculated at the edge of the plate is insensitive to the value of the aspect ratio since it is almost constant for all the aspect ratios studied for the same Rayleigh number. This observation is in good agreement with Boetcher and Sparrow (7) .

Mean Heat transfer
Variations of the mean Nusselt number versus Rayleigh number are computed using Eq. (14). Figures 8 and 9 show the variations of Nusselt number versus Rayleigh number for aspect ratio values ranging from 0.1 to 1.  It can be deduced From Figures 8 and 9 that, for both plates, when Ra increases Nusselt number increases in two stages. The first one is characterized by a low rate of variation due to the dominance of the conduction mode. The second stage is indicated by an important variation of Nusselt number which is enhanced by the onset of convection flow. This transition is performed when Ra exceed a value denoted Ra c . This value is determined graphically from the curve that represents ψ max − ψ min in the space https://www.indjst.org/ between the circular plates versus Rayleigh number as mentioned in Figure 10. The slope changes of (ψ max − ψ min ) emphasizes that the heat transfer within the cavity becomes enhanced by the fluid motion. Similar behavior is observed for Rayleigh-Bénard convection in closed cavities or cavities formed by two infinite (1) plan heated from below where the onset of convection, in the form of rotating fluid cells, starts when Ra number exceed 1708 independently of aspect ratio number. But in our case where the cavity is open, the Ra c that characterize the onset of convection depends on aspect ratio values. Indeed, more the space between the plates is narrowest more Ra c is higher. Figure 11 shows the variations of Ra c in function of A. This graphical representation leads us to model Ra c by a mathematical expression in a power law, the coefficients of the correlations have been computed using the least square method leading to the following expressions: In order to quantify the heat transfer, useful correlations of Nusselt number for the heated plate has been performed in previous works (Table ?? ). In this paper we propose, for the first time in our knowledge, to correlate, the average Nusselt https://www.indjst.org/ number to the Rayleigh number and aspect ratio as a power law for the upper cold plate. The coefficients of the correlations have been computed using the least square method leading to the following expression: The last proposed equation Eq. (16) is plotted in Figure 9 against the numerical results in order to show their accordance. Table 6 shows the aspect ratio and the Rayleigh number ranges where the correlation is valid. The relative error of the correlated data is calculated for every point (A, Ra) with the following expression: Eq. (16) is plotted in Figure 12 against Eq. (18) the correlation given in (18) for the lower plate and the Eq (19) the correlation performed in (5) for upper heated plate facing an adiabatic plan.
Comparing Eq (16) and Eq (18), it can be noted that as Ra number starts to increase . This can be explained by the heated fluid that rises near to the cooled upper circular plate then an important vertical temperature gradient is established in the vicinity of this surface. The same phenomenon was noted by Vafai and Eiteffagh in (4) in an open-ended rectangular cavity with the top and the bottom walls are heated at the same uniform temperature. Figure 12 shows proportionality between − Nu(low) obtained in (18) for the heated plate and − Nu(up) given in (5) . The proportionality coefficient grows with the increase of aspect ratio. In this comparison,  (18) and the correlation − Nu(up) given in (5) .

Conclusion
The free convective heat transfer characteristics, between two circular plates for different values of Rayleigh number and aspect ratios were investigated. A control volume-based finite differences method was used to solve the governing equations numerically. The main conclusions of this study are as follows: • Heat transfer increases with the increase of Rayleigh number and aspect ratio. • The local Nusselt number and its maximum which is located at the edge of the plate increase with the Rayleigh number. • In the case of low Rayleigh numbers or for the low aspect ratios in the case of moderate Rayleigh numbers, the local Nusselt number is constant throughout the surface except at the edge of the plate where it takes a maximum value. • The maximum values of Nusselt number calculated at the edge of the plate is insensitive to the value of the aspect ratio. • For low aspect ratios under moderate to high Rayleigh numbers, the variation of the Local Nusselt number is periodic or pseudo-periodic in the radial direction. The amplitude of these oscillations increases while Rayleigh number increases. • The onset of convection flow is performed at a critical Rayleigh number Ra c that depends on aspect ratio. Eq. (15) is elaborated in order to calculate Ra c .