Convective thermal losses of long-term underground hot water storage
Journal of Thermal Engineering 2024, Vol. 10, Issue 2, pp. 490-502; doi.org/10.18186/thermal.1457058
Abstract
Keywords: Closed Cavern; Heat Transfer; FDM Numerical Solution; Viscous Fluid; Vorticity-Stream Function Formulation; Unsteady Natural Convection
Introduction
Thermal energy storage is playing a vital role in various engineering applications in the past decades and its role is expected to increase in the years to come [1]. Thermal storage systems allow the utilization of waste heat [2], as well as harvesting of renewable solar energy whenever no district heating grid is available or grid capacities are exceeded [3]. They can also balance the power grids through a combined use with heat pumps [4] and along with other storage technologies plays an important role in the smart grid transition of this century [5]. Furthermore, thermal accumulation is crucial for the energy consumption in the building sector [6], especially for buildings with low thermal mass and accumulation capacities [7]. While some recent studies
focus on the incorporation of phase-change materials [7,8], the most mature and applicable technology remains the utilization of sensible thermal energy in water tanks and basins and especially the underground facilities [9] which save urban space and can be incorporated in existing or planned infrastructure. Comparison of state of the art storage technologies acknowledges the sensible thermal energy storage (STES) as one of the cheapest and most common historically proven solutions [10]. A summary of the energy density of the most common STES types shows that the hot water storage is the one with the highest heat capacity per m3 from the sensible heat storage technologies, namely: 60-80 kWh/m3 [11]. Fernandez et al. [12] give a comparison of some of the most popular
*Corresponding author. *E-mail address: mrashevski@imbm.bas.bg This paper was recommended for publication in revised form by Editorin-Chief Ahmet Selim Dalkılıç Published by Yıldız Technical University Press, İstanbul, Turkey Copyright 2021, Yıldız Technical University. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
materials and explain the methodology to assess their practicability in terms of cost and volume. A recent review [3] stresses on the water as best suited storage medium for underground thermal energy storage and its particularly high usability in district heating systems. Consul et al. [13] deal with numerical simulation procedures for pressure-velocity coupling of the Navier-Stokes equations in Boussinesq approximation to compare them with experimental data for a relatively small storage device and a short time period. The interest to correctly understand and model the stratified STES remains a point of discussion as seen in the recent comprehensive review [14] assessing advantages and drawbacks of current analytical and numerical methods for evaluation of thermal energy storage systems. In order to understand better the processes occurring in the hot water storage and the related convective thermal losses, we need to go back to some of the classical works, such as the one by Ostrach [15] - ”Natural convection in enclosures”. The vorticity transport equation is well known and profoundly studied approach to solve the incompressible Navier-Stokes equations, which eliminates the pressure gradient from them. In his works Gresho [16,17] compares the velocity-pressure and the vorticity-stream function formulations and shows that the main challenge of the former one is the determination of accurate vorticity boundary conditions. The most common vorticity boundary conditions are discussed, including Thom’s formula [18] and Wilkes’ one (or sometimes referred as Wilkes-Pearson’s), but also the ones derived by Orszag and Israeli in their comprehensive review [19]. Weinan and Liu [20] give a detailed analysis of the vorticity boundary conditions for FDM, incl. some considerations on the stability issues. A review of some numerical procedures for the vorticity-stream function equations is presented by Tezduyar et al. [21]. A recent work [22] stresses on the importance of the formulation of the vorticity boundary conditions on non-slip walls for the stability of the numerical solution. De Vahl Davis [23] and Saitoh [24] give an accurate benchmark solutions for natural convection in a squared cavity using FDM procedure. Unsteady convection flows for high values of the Rayleigh number are studied using different numerical procedures for short periods (e.g.[13], [25]). Noteworthy is the work of Papanicolau and Belessiotis [26] which aims at providing insight into the behavior of the system at the boundary between laminar and turbulent flow so that the appropriate numerical treatment may be adopted. They use experimental data for the validation of computer codes of two-dimensional buoyant cavity flows [27], and show a good agreement with the model predictions for both the laminar and turbulent regimes. In their study they found out that laminar flows are present up to Ra = 1013, while turbulent regime is expected for Ra ≥ 5x1013. Later on, Hmouda et al. [25] use a different numerical approach to calculate the cooling process by laminar natural convection in cylindrical cavities, once again showing good match with experimental device up to Ra = 3.8x1012. It is easy to see
that with the cooling process, heat losses through the vessel enclosure decrease, effectively reducing the Rayleigh number. Considering this, convection flows up to Ra = 1012 are assumed to remain laminar. In the present work we consider the cooling process of a rectangular two-dimensional domain with time-dependent temperature boundary conditions, describing in a realistic way variable heat losses to the natural environment of an underground sensible thermal energy storage. Flow convection and heat transfer are simulated numerically using FDM in an exemplary design case over 6 months. The numerical procedure is validated with the benchmark solution of De Vahl Davis [23] for natural convection in a square cavity. Solutions of the vorticity and stream function equations are obtained for the Rayleigh numbers Ra = 109 and Ra = 1010 and the results are compared to each other. The study depicts the influence of the Rayleigh number on the creation of vortices and the thermal stratification. Thermal performance of the storage tank (available energy and net efficiency) at the end of the time period and the evolution of the thermocline thickness are presented. Average Nusselt numbers on the horizontal and vertical boundaries of the cavern are calculated. The model gives accurate numerical results with limited computational resources and can be applied for other rectangular shapes, also perfectly applicable for three dimensional facilities with symmetric base.
Formulation Of The Problem
The present work deals with unsteady laminar convection in a closed plane cavern of width d and height H, shown in Figure 1. The cavern is filled with water with an initial temperature Tinit which gradually cools down for a 6 month period without further heat sources. The fluid starts motion that develops in time. The convection is described by the unsteady Navier-Stokes-Fourier equations in Boussinesq approximation, written in dimensionless vorticity-stream function form [28]: (1) (2)
(3) where Pr = ν/κ is the Prandtl number and Ra = gβd3ΔT/ κν is the Rayleigh number. The stream function ψ and vorticity ω are defined by: (4)
Here u and v are x and y velocity components, ν kinematic viscosity, κ = k/ρCP thermal diffusivity, k thermal conductivity, CP specific heat capacity, β thermal expansion coefficient, T dimensionless temperature. The system equations are normalized using scales d, κ/d , ρ, d2/κ and ΔT = Tinit–Tinf for length, velocity, density, time and temperature, correspondingly. Note that the characteristic velocity is caused by the thermal diffusion inside the cavern and the dimensional temperature is measured from the mean temperature Tav = (Tinit+Tinf ) /2. For sufficiently long periods of time the problem studied will turn to a stationary state, which leads to the expectation that at certain moment the temperature in the cavern will fall down to Tav, thus, the deviation from it would not be excessive and will vary between -0.5 and 0.5. The boundary and initial conditions for the stream function and the vorticity are:
(8) The boundary conditions for the temperature are changing at each time step. The temperatures on the walls are calculated using the thermal balance between the walls and the surrounding environment taking into account different insulation thicknesses - d1 at the bottom wall, d2 at the side walls and d3 at the top wall with thermal conductivity kins. The thermal resistance is calculated accordingly as: (9)
Temperature Boundary Conditions The initial condition for the fluid temperature is:
Here R1, R2 and R3 correspond to the thermal resistances for the three different insulation thicknesses, RGR is the mean thermal resistance of surrounding soil layer with thickness dGR and conductivity kGR. We define RW as the sum of conductive thermal resistance of water with thickness dx = x2–x1 and the corresponding average convective thermal resistance of the boundary layer 1/hw. Heat fluxes through the four boundaries are calculated using Fourier equation – three of them facing earth with constant temperature Tinf and the top one facing exterior air with constant temperature Tair. At the top wall, the thermal resistance for exterior boundary air layer, RSE, is accepted in accordance with ISO6946 [29]. On Figure 2 the left side wall composition is shown with the corresponding heat fluxes and temperatures. The boundary temperature on the interior surface of the water storage T1 is found by equalizing the heat fluxes on both
Figure 1. Underground hot water storage at the initial state t=0.
Figure 2. Left side wall composition used to calculate the current boundary temperature T1.
(6) (7) where H̃ =H/d is the dimensionless height of the cavern.
sides of the cavern wall. The unknown temperature T0 is substituted in the system of equations and using the calculated value T2 for each time step, the new interior surface boundary temperature T1 is obtained. The same procedure can be applied to the bottom and top walls. If we introduce a new boundary function fi = RW / (RW+Ri+RGR), we can express the temperature boundary condition for each time step as:
A grid independence check was performed increasing the grid mesh up to 11011 points (at 91x121 nodes). The change in the grid mesh showed only minor variations in the studied variables such as 0.3 % change of the minimum temperature,
0.1. % change of the maximum temperature and 3.6 % change
of the maximum velocity vector after more than 180 days of calculation period, which indicates that the chosen mesh size is adequate and stable. On the other hand, the reduction of the grid points (e.g.31x41 nodes) leads to some instability issues and discrepancies of the results obtained.
Stability and Validation The stability of the whole solution depends strongly on the correct formulation of the vorticity boundary conditions as stated by Borah [22]. For rectangular domains, the classical formulas for the vorticity on the walls calculated with FDM, are summarized by Weinan and Liu [20]. In their work they describe the von Neumann stability condition of the numerical scheme in the interior points which depends on the time and space steps according to:
These formulas define the interior boundary temperatures T1 at each moment t > 0. Numerical Solution Equations (1-3) are solved by explicit FDM, discretizing them in time and space. The Poisson equation of the stream function (3) is solved first using successive over-relaxation (SOR). The boundary conditions for the vorticity are discretized using the Wilkes formula cited by Roache [30], Gresho [17] and Weinan et al. [20] as e.g. at the left wall: (11) The accuracy of this formula for the vorticity at the walls and the solution in the cavern interior points have both an error magnitude of O(h2). As soon as the stream function is calculated and the vorticity boundary conditions are determined, the vorticity diffusion equation is solved. Finally, the temperature equation (2) is solved and the temperature boundary conditions on the four boundaries are updated as explained in the previous section. The algorithm stops when it reaches the final time step. The calculated temperature and velocities are returned to real values in dimensional form and plotted. The calculation time strongly depends on the convergence of the iterations of the stream function. Most of the iterations are performed until accuracy 0.001 is reached. The increase of the accuracy 10 times leads to a significant increase in the calculation time, while some integral characteristics such as the average Nusselt number or the mean storage temperature at the end of the period differ only after the fifth decimal place. The convergence depends on the SOR parameter, being faster when it is closer to one or zero, and slower when it approaches 0.5. Furthermore, the convergence depends on the Rayleigh number. The higher the number, the slower the calculation procedure. On a computer processor Intel Core i7-4700MQ 2.4GHz with 64x OS and with a given accuracy of 0.001 for a mesh grid consisting of 4941 (61x81 points) the calculation times for the total period of 180 days vary between 1 (Ra = 109) and 4 hours (Ra = 1010) with time step of 120 seconds.
(12) Another stability criteria requires restrictions on the flow characteristics and binds the extremum of the velocities to the space time discretization through the following formula: (13) While the first stability criteria is a static one with the thermal diffusivity assumed as constant physical parameter, the second condition requires further dynamic proof throughout the studied period. Figure 3 shows how the function (13) varies over time. Due to the higher velocities present in the case of Ra = 1010, one can observe this case to be closer but still under the threshold value of 2.
Figure 3. Stability criteria C(t) over the entire time period for Ra = 1010 (blue) and Ra = 109 (orange).
Additional check of the calculation procedure is made using the benchmark solutions presented by De Vahl Davis in [23] with the proper boundary conditions. The numerical solution of the present work for Ra = 106 and Pr = 0.71 is compared with the results obtained in his publication and matches them with insignificant margin of error. Some physical characteristics such as the maximum values of the stream function, the velocities and the average Nusselt number (14)
Tinf = 288.15K (15oC), Tair = 273.15K (0oC), H = 8m, d = 6m, and the corresponding thermal properties of water at 363 K: ν = 3.35 x 10−7[m2/s ], κ = 1.66 x 10−7[m2/s ], CP = 4207 [ J/kgK], k = 0.67 [ W/mK], β = 7 x 10−4[K−1]. The bottom boundary insulation thickness is d1 = 0.7m, the side walls have insulation thickness d2 = 0.8m, and the top wall: Table 1. Comparison of some physical characteristics as described by [23] for Ra = 106 and Pr = 0.71 after steady state is reached.
are presented in Table 1. Graphics for the streamlines and the temperature fields of these flows (Figure 4) match well with the ones shown in the benchmark paper [23].
Case Study And Results
Results are plotted in dimensional form for the case studied. The specific parameters used are Tinit = 363.15K(90oC),
(b) (a) Figure 4. Streamlines (4a) and temperature field (4b) for Ra = 106 and Pr = 0.71 after steady state is reached. Streamline contours at -16.79 (0.56) 0 and temperature contours at 0 (0.05) 1.
Figure 5. Time dependent minimum (orange) and maximum (blue) of the dimensionless stream function ψ.
d3=1.2m. The insulation considered is extruded polystyrene foam (XPS) with thermal conductivity kins=0.032[ W/mK]. The dimensionless parameters used for the calculations are Pr = 1.965, corresponding to the given initial temperature of water, and Rayleigh number Ra = 109 to 1010, which correspond to laminar convection regime. As explained previously and shown by the theoretical and experimental studies [26, 27], the flow is expected to remain laminar for values higher than the ones studied by at least two orders of magnitude. Velocity and Temperature Fields The mesh density used to calculate the velocity and temperature fields corresponds to a space step of 10 cm and the equations are calculated with a time step of 2 min. The process starts from an uniform initial temperature, which gradually cools down over a period of 180 days. Once the results are obtained, the values are returned to real values and plotted in the desired frame rate. Extract from all results are shown in Figures 5-11. The maximum and minimum of the stream function fluctuate indicating the strength and number of vortices (Figure 5). The closer their values are – the more, but less intensive vortices are present. The peaks (counterclockwise movement) and the dips (clockwise movement) correspond to the formation of fewer, but bigger vortices. The maximum and minimum oscillate around some mean values except for the beginning of the period when in all cases some significant peaks/dips are observed. Figures 6-11 show the effect of the Rayleigh number on the convection - the higher the value, the more dynamic becomes the fluid motion inside the cavern. The areas of minimum velocity are located in the center of the vortices regardless of their orientation as well as in the bottom tank area (Figures 6-7). The maximum length of the velocity vector Wmax approaches 23.5 cm/min for the higher Rayleigh number in the first third of the studied period. Some values of the maximum velocity at different moments of time are shown in Tables 2-3. The correspondence between the velocity field and the stream lines is well observed comparing Figures 6 and 7 with Figures 8 and 9 accordingly. The negative values of the stream function are marking clockwise vortex orientation,
while the positive ones – counterclockwise. After certain period of time (approx. 29 days at Ra = 109 and 7 days for Ra = 1010) the flow symmetry is lost. Vortices start to build up along the walls and in the upper part of the cavern – counterclockwise (positive) on the left wall and clockwise (negative) at the right wall. Both orientations are present at the top wall. When vortices with the same orientation collide, they get strength and form one bigger vortex. When vortices with different orientation approach and mix up, the movement relaxes and they lose intensity. The correspondence between the thermal and the dynamic pictures is clearly visible - the vortices pull down streams of colder fluid inside warmer field and mix up the water (Figures 10-11). Generally, the temperatures observed in the lower part of the cavern are lower than in the upper half. At the end of the period, the movement starts to calm down and the fluid stratifies vertically along the cavern with more clear vortices at the upper part of the cavern and calmer lower part with less motion. Thermal Performance of the Storage Tables 2 - 3 show the maximum, mean and minimum temperatures inside the cavern and the maximum length of the velocity. With the progress of time, the difference between maximum and minimum temperatures increases reaching 1.9oC for Ra = 109 and 2.1oC for Ra = 1010 after 180 days. The mean temperature of the water is calculated over all mesh-grid points. The difference between the initial water temperature and the mean water temperature in different time steps shows the amount of thermal losses throughout the studied period, taking into account the specific heat capacity of the water and the total volume of the tank. For example, for Ra = 109 the total thermal losses after 180 days is 441.4 kWh, while for Ra = 1010 for the same period is 559.0 kWh or 21% higher than the previous case. This shows the importance of the fluid convective behavior for the seasonal efficiency of a well-insulated hot water storage. The mean water temperature is important to calculate the available thermal energy of the water storage according to the following formulas:
Table 2. Temperature (oC) and maximum velocity (cm / min) values at different days for Ra = 109
Table 3. Temperature (oC) and maximum velocity (cm / min) values at different days for Ra = 1010
Figure 6. Velocity fields for Ra = 109 at different time stops. Figure 7. Velocity fields for Ra = 1010 at different time stops.
(15) where TH is the delivery temperature for the heating system and V is the total volume of the water storage. The lower the delivery temperature for the heating system (e.g. hydronic radiant heating systems at TH = 35o C), the bigger the temperature difference, the higher the amount of available thermal energy in the storage facility at the end
of the period. Assuming low temperature radiant hydronic heating, the amount of energy in the storage after 180 days is 60.1 kWh/m3 (in accordance with the numbers given by Schmidt [11] - between 60 and 80 kWh/m3), which can be increased further if a heat pump is added to the system design. The efficiency of the storage is calculated with the formula ηD = Qfin /Qinit (see [14]). Depending on TH (e.g. between 35oC and 60oC) the values for ηD vary between 51 and 97%. Even though the design of the tank and the
Figure 9. Streamlines for Ra = 1010 at different time stops.
operation conditions are different, these results can be formally compared with the underground seasonal pit storage in Dronninglund [31] showing efficiency levels between 90 and 96% and maximum temperatures between 84 and 89oC. The convection losses are lower for the lower value of the Rayleigh number, thus, the gap between maximum and minimum temperatures increases more slowly and the mean temperature of the water remains higher for longer period of time.
Another indicator for the thermal stratification of the storage tank is the thermocline thickness. The thermocline thickness is the region of the tank where strong change of the temperature in vertical direction is observed. According to Bahnfleth and Musser [32], it is the region, where 0.1 ≤ Θ ≤ 0.9 and Θ is defined as: (16)
Figure12 shows the thermocline profiles in the center of the cavern along its height for Ra = 109 and Ra = 1010 in time intervals of two months. The more stratified is the water, the smaller is the thermocline thickness, as shown by Haller et al. [33]. In our case, however, no clear cold water zone is visible, because there is no cold water inlet from below. We can easily recognize that with the time, the water in the storage tends to mix up and the thermocline
thickness increases. For example, for Ra = 1010 it increases approximately from 2 to 3m height between the second (Θ ≤ 0.80) and the sixth month (Θ ≤ 0.85). Nevertheless, the hot water zone is clearly visible, decreases very slowly and remains in more than the half of the tank height even after 6 months period. The fluid temperature in the very top of the cavern decreases almost linearly due to direct thermal losses through the upper boundary. For the higher Rayleigh
Figure 10.Temperature fields for Ra=109 at different time Figure 11.Temperature fields for Ra=1010 at different time stops stops.
Figure. 12. Θ in the middle of the water tank (x = 0.5) for days 180 (blue), 120 (red) and 60 (yellow). number the water mixing diminishes the thickness of this boundary layer. Nusselt Number Unlike the benchmark solution, where the heat flux is directed from the left to the right side wall only, in our case it is going through all the four walls. Since the average Nusselt number presents the ratio of the average heat flux between the convective and conductive cases, its values should be calculated after solving the pure conductive thermal problem for solid body cooling (Figure 13). For example, on the left/right side wall the Nusselt number is: Figure. 13. Temperature field TCOND for solid body cooling.
Time-dependent graphics are presented for the values of the Nusselt number for the bottom, left/right side and top walls of the cavern (Figure 14). The values of the Nusselt
Figure 14. Nusselt number for Ra = 10 (a) and Ra = 10 (b) for the top (1), left/right side (2) and bottom (3) walls. 9
number on the left and right side walls are, as expected, equal to each other, because the boundary conditions are the same in any moment of time. The calculation of this integral characteristic is an additional check of the stability of the numerical scheme. The curves for the vertical and upper walls are fluctuating stronger than for the bottom one, because of the intensity of the vortices observed in proximity to those walls. The higher the Rayleigh number, the stronger the mixing process in the cavern, the smaller the convective heat flux near the walls. This results in slightly steeper curves with the progress of time. For both Rayleigh numbers at the beginning of the studied period the Nusselt number increases on the top and side walls, marking the intensity of the convective heat transfer. While the Nusselt number on those walls for Ra = 109 reaches its maximum at the end of the studied period, for Ra = 1010 it passes through its maximum and declines. This declination is occurring earlier for the higher Rayleigh number, due to the faster cavern cooling and decreased temperature difference at the walls leading to reduced convective heat transfer. In the bottom part of the cavern the convective heat flux is lower than the conductive one, thus, the Nusselt number is less than one. This effect is due to the cool fluid layer with less fluid motion appearing in the lower part of the cavern as visible from the temperature and streamline profiles (Figures 8-11).
Conclusions
Flow convection and heat transfer in a closed rectangular two-dimensional cavern are studied. Numerical solutions for the stream function equation, the vorticity transport equation and the temperature equation are obtained using FDM. The temperature equation is solved with time-dependent temperature boundary conditions on the cavern walls formulated through Fourier’s law. The numerical procedure satisfies two stability criteria and is validated by a comparison with other works. The procedure showed good matching with the benchmark results for natural convection in a quadratic cavity obtained by previous authors. The convective processes are analyzed through a case study of an exemplary underground hot water storage with given geometrical and thermal properties. The storage tank is investigated for a six month period. Solutions are obtained for two values of the Rayleigh number and the corresponding graphics for the velocity, stream function and temperature fields are presented. The effect of the Rayleigh number on the convective vortex formation, absolute maximum length of the velocity vector, temperature distribution and thermal losses is shown through a comparison between the results for Ra = 109 and Ra = 1010. The values of the stream function are representative for the number of vortices and their intensity: when the maximum and minimum values are closer the vortices are more in number but less
in intensity. If the extreme values of the stream function are further separated, the vortices are fewer but more intensive. Collision of vortices with the same orientation forms bigger and more intensive ones, while collision of vortices with opposite orientation reduces their size and intensity. The cooling of the cavern over long periods of time is graphically and numerically presented through the extreme and mean values of the temperature field. Starting from an unbalanced state with an uniform initial temperature, the thermal stratification becomes obvious with the progress of time. The variation of the thermocline thickness is shown to be comparatively small, being stronger for the higher Rayleigh number due to the stronger mixing processes. Time-dependent values of the average Nusselt number on the bottom, top and side walls of the cavern are determined and presented graphically. The declination of the Nusselt number curves at the end of the period occurs earlier for the higher Rayleigh number, due to the stronger mixing processes. As a result, the convective heat transfer through the boundaries slows down with time, due to the decreased temperature difference at the walls. The small thermal losses and the available thermal energy at the end of the period indicate good potential for long-term energy storage.
Nomenclature
C CP d d1, d2, d3 dGR g H h hW k kins kGR Nu Pr qconv, qcond Qinit, Qfin R1, R2, R3 RGR RW Ra T T0 Tair , Tinf Tav
Second stability criteria [-] Specific heat [J/kgK] Width of the cavern [m] Wall insulation thicknesses [m] Thicknesses of the adjacent ground layer[m] Gravitational acceleration [m/s2] Height of the cavern [m] Discretization space step [-] Heat transfer coefficient [W/m2K] Fluid thermal conductivity[W/mK] Conductivity of the thermal insulation [W/mK] Conductivity of the adjacent ground layer [W/mK] Nusselt number [-] Prandtl number [-] Dimensionless heat fluxes through the walls [-] Available thermal energy [J] Wall thermal resistances [m2K/W] Ground thermal resistance [m2K/W] Thermal resistance of water [m2K/W] Rayleigh number [-] t Time [s] Fluid temperature field[K] Reference temperature [K] Exterior boundary temperatures [K] Mean temperature [K]
Temperature field for solid body cooling[K] Initial fluid temperature [K] Heating system temperature [K] Storage fluid temperatures [K] Velocity components [m/s] Maximum length of the velocity vector [m/s] Cartesian coordinates [m]
Thermal expansion coefficient [1/K] Laplacian operator Thermal efficiency [-] Thermocline temperature ratio [-] Thermal diffusivity [m2/s] Kinematic viscosity [m2/s] Stream function [m2/s] Vorticity [1/s]
Acknowledgements
This work was supported by the Bulgarian Ministry of Education and Science under the National Research Programme “Young scientists and postdoctoral students -2” approved by DCM 206 / 07.04.2022.
Data Availability Statement
The authors confirm that the data that supports the findings of this study are available within the article. Raw data that support the finding of this study are available from the corresponding author, upon reasonable request.
Conflict Of Interest
The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Ethics
There are no ethical issues with the publication of this manuscript.
References
- Heat Transf 1972;8:161−227. [CrossRef]
- Chavan S, Rudrapati R, Manickam S. A compre- [16] Gresho PM. Some current CFD issues relevant to the hensive review on current advances of thermal incompressible Navier-Stokes equations. Comp Math energy storage and its applications. Alexandria Eng Appl Mech Engineer 1991;87:201−252. [CrossRef] J 2022;61:5455−5463. [CrossRef] [17] Gresho PM. Some interesting issues in incom-
- Miro L, Gasia J, Cabeza LF. Thermal energy storage pressible fluid dynamics, both in the continuum (TES) for industrial waste heat (IWH) recovery: A and in numerical simulation. Adv Appl Mech review. Appl Energy 2016;179:284−301. [CrossRef] 1991;28:45−140. [CrossRef] 502 J Ther Eng, Vol. 10, No. 2, pp. 490−502, March, 2024
- Thom A. The flow past circular cylinders at low [26] Papanicolaou E, Belessiotis V. Transient natural con- speeds. Proc Royal Soc London 1933;A 141:651−666. vection in a cylindrical enclosure at high Rayleigh [CrossRef] numbers. Int J Heat Mass Transf 2002;45:1425−1444.
- Orszag SA, Israeli M. Numerical simulation of vis- [CrossRef] cous incompressible flows. Annu Rev Fluid Mech [27] Cheesewright R, King K, Ziai S. Experimental data 1974;6:281−318. [CrossRef] for the validation of computer codes for the predic-
- Weinan E, Liu J. Vorticity boundary condition and tion of two-dimensional buoyant cavity flows. Am related issues for finite difference schemes. J Comp Soc Mech Eng Heat Transf Div 1986;60:75−81. Phys 1996;124:368−382. [CrossRef] [28] Jaluria Y. Natural Convection. Heat and Mass
- Tezduyar T, Liou J, Ganjoo D, Behr M. Solution tech- Transfer. Oxford: Pergamon Press; 1980. niques for the vorticity-stream function formulation [29] ISO 6946:2017. Building components and building of two-dimensional unsteady incompressible flows. elements - thermal resistance and thermal transmit- Int J Num Meth Fluids 1990;20:515−539. [CrossRef] tance - calculation methods Standard. International
- Borah A. Computational study of stream func- Organization for Standardization, Geneva, CH. 2017. tion-vorticity formulation of incompressible flow [30] Roache P. Computational Fluid Dynamics. and heat transfer problems. Appl Mech Math Albuquerque: Hermosa Publishers; 1972. 2011;52-54:511−516. [CrossRef] [31] Winterscheid C, Schmidt T. Dronninglund district
- De Vahl Davis G. Natural convection of air in a heating monitoring data evaluation for the years square cavity. A bench mark numerical solution. Int 2015-2017. Stuttgart: Solites; 2019. J Num Meth Fluids 1983;3:249−264. [CrossRef] [32] Bahnfleth W, Musser A. Thermal performance of
- Saitoh T. High-accuracy bench mark solutions to a full-scale stratified chilled-water thermal storage natural convection in a square cavity. Comp Mech tank. ASHRAE Transf 1998;104:377−388. 1989;4:417−427. [CrossRef] [33] Haller M, Cruickshank C, Streicher W, Harrison S,
- Hmouda I, Rodriguez I, Bouden C, Oliva A. Andersen E, Furbo S. Methods to determine stratifi- Unsteady natural convection cooling of a water stor- cation efficiency of thermal energy storage processes age tank with an internal gas flue. Int J Therm Sci - review and theoretical comparison. Solar Energy 2010;49:36−47. [CrossRef] 2009;83:1847−1860. [CrossRef]
Share and Cite
RASHEVSKI, M.; SLAVTCHEV, S. Convective thermal losses of long-term underground hot water storage. Journal of Thermal Engineering 2024, Vol. 10, pp. 490-502. https://doi.org/10.18186/thermal.1457058