Transient Natural Convection in a Thermally Insulated Annular Cylinder Exposed to a High Temperature from the Inner Radius

: Extensive numerical analysis was performed for the unsteady state, natural convection in the annular cylinders. The cylinder’s boundaries were thermally insulated, except the inner surface. The ﬂuid (water) in the cylinder initially was assumed at a cold temperature while the inner surface was subjected to a high temperature. The time history for the heat transfer by di ﬀ usion and advection was studied. The time needed for fully charging the storage tank and rate of heat transfer was calculated. The predicted results were compared with the pure heat di ﬀ usion process and with a steady-state convection system. Therefore, CFD simulations were performed for natural convection in the storage tank. The main objective of this study was to establish correlations for the rate of heat transfer as a function of time and other controlling parameters. The correlation is needed in designing a thermal energy storage system for domestic and industrial heating processes. One of the drawbacks of the conventional thermal storage systems is the slow charging and discharging, where the heat transfer is mainly di ﬀ usion dominated. To overcome such a problem, a system was designed based on the natural convective heat transfer mechanism. Therefore, the heat transfer and ﬂuid ﬂow in a cylindrical storage tank were simulated for a range of Rayleigh numbers (10 4 to 10 8 ) and radius ratio. It was found that a convection-operated storage tank reduces the thermal charging process time drastically compared with the thermally di ﬀ usion charging process. The rate of reduction in the charging time mainly depends on the rate of heating and geometric parameter of the tank. To the best of the authors’ knowledge, the work is novel.


Introduction
Thermal storage systems are an essential part of domestic and commercial water heating systems, especially solar heating systems, for water heating and power generation. It is crucial to understand the flow and heat transfer in sensible and latent heat thermal storage systems in designing and optimizing the storage system. One of the critical parameters in the design of the system is estimating the charging and discharging times. The steady and unsteady states of natural convection in enclosures have been extensively studied by Hussein et al. [1].
However, very limited data are available on the natural convection in annular geometries, specifically unsteady state convection with thermally insulated boundaries. The operation of thermal storage systems is unsteady due to charging and discharging processes. From the physics of heat transfer, in the very early thermal charging process, the expected mechanism of heat transfer is the could cause instability of the equilibrium state of magnetic nanofluid and influence gravity convection. The authors considered Prandtl numbers from 0.7 to 700 and Rayleigh numbers of up to 10 6 . The authors found that the magnetic permeability of the magnetic fluid filling the enclosure was an important factor. They studied four modes of fluid flow analytically, and two of them were demonstrated using the numerical simulations. The simulation showed that the outer uniform magnetic field turned to nonuniform in the magnetic fluid heated from below. So, there was a weak convective flow, even for a very small magnetic Rayleigh number.
Li et al. [8] studied the effect of different thermal storage tank structures on temperature stratification and thermal efficiency during the charging process. The authors compared temperature stratification in six kinds of tanks. Also, the thermal charging efficiency of the tanks was compared. The optimal tank shape for charging was developed. The authors also performed CFD simulation and validated their results experimentally and found the relative discrepancy between the CFD model predictions and experimental data was about 5%.
Medebber and Retiel [9] examined the effects of buoyancy and height ratio on the natural convection in a vertical cylinder partially annular. The results showed that the flow field is characterized by simple recirculating zone turning in the clockwise direction. The results indicated that the variation in the Ra number and height ratio have a major influence on the flow structure and isotherm patterns. They noticed two flow regimes for lower Ra number values and the dominance of conduction heat transfer. At higher Ra number values, the heat transfer rate was increased and was dominated by convection mode. Moreover, the results showed that the average Nusselt number increased with the increasing of the Ra number and height ratio, which is consistent with the predications of previously mentioned authors.
Moran and Katz [10] numerically investigated a natural convection heat transfer in an annular cylindrical enclosure with nonuniform boundary conditions using the finite element method. The inner circumferential temperature distribution was assumed to be varied along the inner hot cylinder surface. The outer surface of the cylinder was kept at a constant cold temperature. The simulation was carried on for Ra = 2 × 10 5 to 1.1 × 10 7 . The authors analyzed the temperature and velocity by comparing the results to the case of uniform boundary conditions. A correction factor was proposed to enable the use of the classical empirical correlations for the uniform boundary conditions. The results showed that the overall heat transfer rate for nonuniform linear temperature distribution was up to 7% higher compared to the cases of uniform temperature equal to the average nonuniform temperature. However, a 7% increase in the rate of heat transfer for an idealized system is not that high for practical applications.
Moldovan et al. [11] experimentally investigated the natural convection flow in an annularly heated vertical cylinder. The work was focused on measuring the flow temperatures and velocities in a crystal grow reactor equivalent setup. The cylinder geometry was scaled to imitate the geometry of existing crystal growth reactors. The Rayleigh number varied from 10 7 to 10 8 were investigated. The PIV was used to measure and visualize the flow velocity. The experimental measurements showed that axial heat transfer occurred mainly by advection, and the radial heat transfer occurred primarily by conduction. Also, they mentioned that 3D cells generated by the temperature gradient improved the heat transfer between the cold and hot sections of the setup.
Bai et al. [12] studied thermal stratification in a cylindrical tank due to heat losses in standby mode. The experimental and numerical simulations were performed. The Nusselt number was determined, and a new one-dimensional model for the cooling process was proposed and experimentally validated. The results showed that the aspect ratios, H/D, higher than three had little influence on the thermal stratification.
Goodrich et al. [13] studied natural convection heat transfer and boundary layer transition for vertical heated cylinders. The effect of curvature on the rate of heat transfer and regime transition was experimentally investigated. The authors applied PIV to visualize the flow field and measured the local surface temperature from five heater sizes. The authors provided a continuous, local, and curvature dependent Nusselt number correlation for the laminar, transition, and turbulent regimes for natural convection from vertically heated cylinders.
The thermal storage systems may be coupled with PV-T panels to utilize the waste heat from PV cooling and heat the domestic water. This kind of system, with tube plate PV/T filled with iron, was studied by Huo at al. [14]. The use of iron-filling allows increasing the electrical and thermal efficiency of the PV-T system.
In conclusion, the mentioned works considered different conditions and scenarios for steady and unsteady state heat convection in enclosures. However, to our knowledge, there has been no work done on the annular cylinder with all surfaces thermally insulated except the heated inner surface. Furthermore, no correlation of the rate of heat transfer for such a system has been reported. Hence, it is essential to investigate such a system to estimate the charging and discharging times and rate of heat transfer. Figure 1 shows the sketch of the problem with coordinate and boundary conditions. The radius aspect ratio R in = r i /r o was set to 0.05, 0.1, and 0.2. The height aspect ratio L* = L/r o was set to one and two. All the boundaries of the enclosure were assumed thermally insulated, except the inner surface of the annular cylinder, which was suddenly heated to a constant temperature while the water in the encloser was at a cold temperature. The water Prandtl (Pr) number was set to 4.0 (water at a moderate temperature). A range of the rate of heating was explored for Ra = 10 4 to Ra = 10 8 . The thermophysical properties of the working fluid were assumed constant, except for the change of the density in the buoyancy term, for which the Boussinesq approximation was employed.

Problem Definition and Modelling Details
Energies 2020, 13, 1291 4 of 13 studied by Huo at al. [14]. The use of iron-filling allows increasing the electrical and thermal efficiency of the PV-T system.
In conclusion, the mentioned works considered different conditions and scenarios for steady and unsteady state heat convection in enclosures. However, to our knowledge, there has been no work done on the annular cylinder with all surfaces thermally insulated except the heated inner surface. Furthermore, no correlation of the rate of heat transfer for such a system has been reported. Hence, it is essential to investigate such a system to estimate the charging and discharging times and rate of heat transfer. Figure 1 shows the sketch of the problem with coordinate and boundary conditions. The radius aspect ratio Rin = ri/ro was set to 0.05, 0.1, and 0.2. The height aspect ratio L* = L/ro was set to one and two. All the boundaries of the enclosure were assumed thermally insulated, except the inner surface of the annular cylinder, which was suddenly heated to a constant temperature while the water in the encloser was at a cold temperature. The water Prandtl (Pr) number was set to 4.0 (water at a moderate temperature). A range of the rate of heating was explored for Ra = 10 4 to Ra = 10 8 . The thermophysical properties of the working fluid were assumed constant, except for the change of the density in the buoyancy term, for which the Boussinesq approximation was employed. Two-dimensional axisymmetric Navier-Stokes equations were solved using the second-order accurate finite-volume method. The governing equations can be written in a nondimensional form.

Problem Definition and Modelling Details
The continuity equation: The momentum equation: R-component: Two-dimensional axisymmetric Navier-Stokes equations were solved using the second-order accurate finite-volume method. The governing equations can be written in a nondimensional form.
The continuity equation: The momentum equation: R-component: Z-component: The energy equation: The outer radius of the cylinder was used as a reference length scale. The nondimensional terms velocity, thermal diffusivity, and outer radius, respectively.
The local Nusselt number was based on the heat balance: where h, k, and R in represent the convective heat transfer coefficient, thermal conductivity of the fluid, and dimensional inner radius of the cylinder R in = (r i /r o ), respectively. The equilibrium amount of heat stored in the enclosure was: Nondimensionalized: where A = L/r o . Also, the equilibrium amount of heat can be calculated as: In the nondimensionalized form: The nondimensional heat stored in the tank as a function of time was: The finite volume method with implicit time discretization scheme was used for solving governing Equations (1)-(4). The computational procedure applied allowed us to obtain the R-velocity and Z-velocity distribution within the tank, as well as temperature distribution within the tank. Those parameters were used to calculate the nondimensional heat stored in the tank as a function of time.
For a fully charged system, Equation (10) should be equal to Equation (8). In the modeling, we used Equation (10) in calculating the amount of energy stored in the system at any time. At the equilibrium state, the prediction of Equation (10) produced the same results of Equation (8). This is an indication that the model correctly predicts the analytical results (Equation (8)). For instance, the dimensionless amount of energy stored for an annular cylinder of R in = 0.1 and A = 1 was 0.99 (Equation (8)), and the model's prediction was 0.98994 (Equation (10)).
The numerical predictions results were carefully tested for grid size and time step independence. The grids were clustered near the boundaries. It was found that for Ra up to 10 6 , the grid sizes of 61 × 61 produced independent results for an aspect ratio of one. Figure 2 shows the grid dependency results for Ra = 10 8 using 81 × 81 and 121 × 121 nonuniform grids. An insignificant difference between the predictions of those grid sizes was noticed. To ensure that the predicted results were grid-independent, 61 × 61 nonuniform grids were used for Ra of 10 4 , 10 5 , and 10 6 . For Ra = 10 7 and Ra = 10 8 , we used 81 × 81 and 121 × 121, respectively. Also, the time step used in updating the results was tested. For instance, for Ra = 10 8 , the time step was set to 1 × 10 −6 . Also, for height aspect ratios more than one, the number of grids increased proportionally along the cylinder height, The average Nusselt number Nu av was calculated by averaging the local Nusselt number determined based on Equation (6) at the location of r = r in .
Energies 2020, 13, 1291 6 of 13 an indication that the model correctly predicts the analytical results (Equation (8)). For instance, the dimensionless amount of energy stored for an annular cylinder of Rin = 0.1 and A = 1 was 0.99 (Equation (8)), and the model's prediction was 0.98994 (Equation (10)). The numerical predictions results were carefully tested for grid size and time step independence. The grids were clustered near the boundaries. It was found that for Ra up to 10 6 , the grid sizes of 61 x 61 produced independent results for an aspect ratio of one. Figure 2 shows the grid dependency results for Ra = 10 8 using 81 × 81 and 121 × 121 nonuniform grids. An insignificant difference between the predictions of those grid sizes was noticed. To ensure that the predicted results were gridindependent, 61 × 61 nonuniform grids were used for Ra of 10 4 , 10 5 , and 10 6 . For Ra = 10 7 and Ra = 10 8 , we used 81 × 81 and 121 × 121, respectively. Also, the time step used in updating the results was tested. For instance, for Ra = 10 8 , the time step was set to 1 × 10 −6 . Also, for height aspect ratios more than one, the number of grids increased proportionally along the cylinder height, The average Nusselt number Nuav was calculated by averaging the local Nusselt number determined based on Equation (6) at the location of r = rin. Moreover, the developed model was tested for an annular cylinder heated from the inner surface and cooled from the outer surface with other boundaries that were assumed adiabatic. The model predictions agree well with the analytical solution.
We can claim with confidence that the predicted results are numerically accurate and reliable.

Results
In most domestic applications, a water temperature of about 45 °C is used. The Prandtl number of water for the mentioned temperature is about 4.0 [15]. It is well established in the open literature that the rate of heat transfer (Nusselt number) for natural convection in enclosures for fluids of Pr > 1 is not a strong function of Pr. Therefore, keeping Pr at 4.0, the results of the current work should be applicable to the working temperature range. Moreover, the developed model was tested for an annular cylinder heated from the inner surface and cooled from the outer surface with other boundaries that were assumed adiabatic. The model predictions agree well with the analytical solution.
We can claim with confidence that the predicted results are numerically accurate and reliable.

Results
In most domestic applications, a water temperature of about 45 • C is used. The Prandtl number of water for the mentioned temperature is about 4.0 [15]. It is well established in the open literature that the rate of heat transfer (Nusselt number) for natural convection in enclosures for fluids of Pr > 1 is not a strong function of Pr. Therefore, keeping Pr at 4.0, the results of the current work should be applicable to the working temperature range.
The typical scenario of fluid flow and heat transfer evolutions are shown in Figures 3 and 4 for streamlines and isotherms, respectively, at different time frames. Figures 3 and 4 are for Ra = 10 6 and for r i * = 0.1. At the early stage, τ = 0.01, shown in Figures 3a and 4a, a weak flow circulation took place inside the enclosure. At this stage, the heat transfer was mainly by diffusion. As time increased, the strength of the flow circulation increased (Figure 3b,c.) However, as the time further precedes the strength of the recirculation starts decaying (Figure 3d). Figure 4 shows isotherm (temperature) evolution at different time frames. At early stages ( Figure 4a) the thermal plume only evidence at the upper part of the cavity. The convection dominated flow was clear in Figure 4b,c. Also, the flow stratification was evident within those time frames. At a dimensionless time of about 0.5 (Figure 4d), the temperature of the fluid tended toward equilibrium. In other words, the system was fully charged, and the temperature was equal to the hot boundary temperature. The typical scenario of fluid flow and heat transfer evolutions are shown in Figures 3 and 4 for streamlines and isotherms, respectively, at different time frames. Figures 3 and 4 are for Ra = 10 6 and for ri* = 0.1. At the early stage, τ = 0.01, shown in Figures 3a and 4a, a weak flow circulation took place inside the enclosure. At this stage, the heat transfer was mainly by diffusion. As time increased, the strength of the flow circulation increased (Figure 3b,c.) However, as the time further precedes the strength of the recirculation starts decaying (Figure 3d). The convection dominated flow was clear in Figure 4b,c. Also, the flow stratification was evident within those time frames. At a time of about 0.5, the temperature of the fluid tended toward equilibrium. In other words, the system was fully charged, and the temperature was equal to the hot boundary temperature.      It is clear that as the rate of heating increased, i.e., Ra increased, the fluid in the enclosure reached the thermal equilibrium condition faster compared with a pure diffusion process. The rate of heat transfer (Nusselt number) can be correlated as (Table 1)   It is clear that as the rate of heating increased, i.e., Ra increased, the fluid in the enclosure reached the thermal equilibrium condition faster compared with a pure diffusion process. The rate of heat transfer (Nusselt number) can be correlated as (Table 1)  It is clear that as the rate of heating increased, i.e., Ra increased, the fluid in the enclosure reached the thermal equilibrium condition faster compared with a pure diffusion process. The rate of heat transfer (Nusselt number) can be correlated as (Table 1):  5 1.366 10 6 2.266 10 7 3.881 10 8 6.753 The coefficient can be correlated as a function of Ra as: For Ra = 10 4 : For Ra = 10 5 , the effect of the radius ratio on Nu can be correlated as: For Ra = 10 6 : For Ra = 10 7 : For Ra = 10 8 :

Energy Storage
The equilibrium thermodynamic helps us to calculate the amount of energy stored in an enclosure, which is given in dimensional and nondimensional forms in Equations (7) and (8), respectively.
However, the equilibrium thermodynamics does not help us in estimating the time of the process. On the other hand, the heat transfer (nonequilibrium thermodynamics) helps us in estimating the time needed for charging the storage system. Figure 8 shows a typical charging process for an annular tank with a radial aspect ratio of 0.1 and the height aspect ratio of one, for different heating processes (Rayleigh numbers). It is very clear that as the heating rate increased, the charging time decreased due to the enhanced buoyancy force. We estimated that the charging time needed for the heat capacity of the storage system reached 99% of the total heat capacity.  Table 2 summarizes the computationally predicted results for a range of investigated Ra and the aspect ratio of one. The convection results are compared with a pure diffusion process, where the flow assumed to be stagnant. For instance, even for low Ra (Ra = 10 6 ), the time needed to charge the enclosure was about 40% of the time needed by the diffusion process.
Some limitations of the model exist. For example, in case of solar air heating, the air temperature decreases along the pipe height L and air inlet temperature changes with time. Therefore, the coupled boundary conditions between the air flowing in the pipe and tank domain should be applied to the model in more detail to conjugate the heat transfer problem [16].

Conclusions
Computational fluid dynamics (CFD) was used to simulate unsteady state natural convection in annular cylinders for a range of Ra up to 10 8 and for radius aspect ratios of 0.05, 0.1, and 0.2. Streamlines, isotherms, and rate of heat transfer (Nusslt number) as a function of time and Rayleigh number were analyzed and reported. The control volume method was employed for discretizing two-dimensional and axisymmetric natural convection problem. The explicit method was used for  Table 2 summarizes the computationally predicted results for a range of investigated Ra and the aspect ratio of one. The convection results are compared with a pure diffusion process, where the flow assumed to be stagnant. For instance, even for low Ra (Ra = 10 6 ), the time needed to charge the enclosure was about 40% of the time needed by the diffusion process. The above data can be correlated as: Some limitations of the model exist. For example, in case of solar air heating, the air temperature decreases along the pipe height L and air inlet temperature changes with time. Therefore, the coupled boundary conditions between the air flowing in the pipe and tank domain should be applied to the model in more detail to conjugate the heat transfer problem [16].

Conclusions
Computational fluid dynamics (CFD) was used to simulate unsteady state natural convection in annular cylinders for a range of Ra up to 10 8 and for radius aspect ratios of 0.05, 0.1, and 0.2. Streamlines, isotherms, and rate of heat transfer (Nusslt number) as a function of time and Rayleigh number were analyzed and reported. The control volume method was employed for discretizing two-dimensional and axisymmetric natural convection problem. The explicit method was used for discretizing time dependent terms of governing equations (continuity, momentum, and energy). The physics of fluid flow and heat transfer in annular heat storage tanks were explained. A few correlations were developed for the rate of heat transfer, which is crucial for designing and operating the thermal storage tanks. The correlations were valid for numbers from 10 4 to 10 8 and allowed us to calculate the Nusselt number changes in time of tank charging. Furthermore, the charging period was analyzed and correlated as a function of Ra number. The work has importance in designing hot water storage tanks for solar energy and other domestic or industrial applications and allows the prediction of the charging time of the tank.