Convective thermal losses of long-term underground hot water storage

A case of underground long-term hot water storage is investigated numerically. The study is based on the unsteady two-dimensional Navier-Stokes equations in Boussinesq approximation applied to a closed cavern with time-dependent temperature boundary conditions on the walls. The problem formulated in a vorticity-stream function statement is solved by finite difference method (FDM) for high values of the Rayleigh number and for the Prandtl number of water. Streamlines, velocity and temperature fields are presented graphically for given moments of time. The evolution of the thermocline thickness in the mid-section of the cavern is slow and illustrates that the hot water zone occupies more than the half of the cavern even after 6 months period. The Nusselt number on the walls shows that the convective thermal losses are small and after certain period of time tend to decrease due to the diminished temperature difference at the walls. The influence of the fluid convection on the thermal losses is evaluated quantitatively, showing high seasonal thermal efficiency of the insulated hot water storage. Cite this article as: Rashevski M, Slavtchev S. Convective thermal losses of long-term under-ground hot water storage. J Ther Eng 2024;10(2):490−502.


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 m 3 from the sensible heat storage technologies, namely: 60-80 kWh/m 3 [11].Fernandez et al. [12] give a comparison of some of the most popular 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 = 10 13 , while turbulent regime is expected for Ra ≥ 5x10 13 .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.8x10 12 .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 = 10 12 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 = 10 9 and Ra = 10 10 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 T init 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]: (3) where Pr = ν/κ is the Prandtl number and Ra = gβd 3 Δ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/ρC P thermal diffusivity, k thermal conductivity, C P specific heat capacity, β thermal expansion coefficient, T dimensionless temperature.The system equations are normalized using scales d, κ/d , ρ, d 2 /κ and ΔT = T init -T inf 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 T av = (T init +T inf ) /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 T av , 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: (5) (6) (7) where H ˜ =H/d is the dimensionless height of the cavern.

Temperature Boundary Conditions
The initial condition for the fluid temperature is: 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 -d 1 at the bottom wall, d 2 at the side walls and d 3 at the top wall with thermal conductivity k ins .The thermal resistance is calculated accordingly as: (9) Here R 1 , R 2 and R 3 correspond to the thermal resistances for the three different insulation thicknesses, R GR is the mean thermal resistance of surrounding soil layer with thickness d GR and conductivity k GR .We define R W as the sum of conductive thermal resistance of water with thickness dx = x 2 -x 1 and the corresponding average convective thermal resistance of the boundary layer 1/h w .Heat fluxes through the four boundaries are calculated using Fourier equation -three of them facing earth with constant temperature T inf and the top one facing exterior air with constant temperature T air .At the top wall, the thermal resistance for exterior boundary air layer, R SE , is accepted in accordance with ISO-6946 [29].
On  sides of the cavern wall.The unknown temperature T 0 is substituted in the system of equations and using the calculated value T 2 for each time step, the new interior surface boundary temperature T 1 is obtained.The same procedure can be applied to the bottom and top walls.
If we introduce a new boundary function f i = R W / (R W +R i +R GR ), we can express the temperature boundary condition for each time step as: (10) These formulas define the interior boundary temperatures T 1 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(h 2 ).
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 = 10 9 ) and 4 hours (Ra = 10 10 ) with time step of 120 seconds.
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: (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 = 10 10 , one can observe this case to be closer but still under the threshold value of 2. 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 = 10 6 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) 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].Table 1.Comparison of some physical characteristics as described by [23] for Ra = 10 6 and Pr = 0.71 after steady state is reached.The dimensionless parameters used for the calculations are Pr = 1.965, corresponding to the given initial temperature of water, and Rayleigh number Ra = 10 9 to 10 10 , 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.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 W max 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 = 10 9 and 7 days for Ra = 10 10 ) the flow symmetry is lost.Vortices start to build up along the walls and in the upper part of the caverncounterclockwise (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.9 o C for Ra = 10 9 and 2.1 o C for Ra = 10 10 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 = 10 9 the total thermal losses after 180 days is 441.4 kWh, while for Ra = 10 10 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: where T H 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 T H = 35 o 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/m 3 (in accordance with the numbers given by Schmidt [11] -between 60 and 80 kWh/m 3 ), 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 = Q fin /Q init (see [14]).Depending on T H (e.g. between 35 o C and 60 o C) the values for η D vary between 51 and 97%.Even though the design of the tank and the   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 89 o C. 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:   Figure12 shows the thermocline profiles in the center of the cavern along its height for Ra = 10 9 and Ra = 10 10 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 = 10 10 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 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: (17) 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    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 = 10 9 reaches its maximum at the end of the studied period, for Ra = 10 10 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 = 10 9 and Ra = 10 10 .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.

Figure 2 Figure 1 .
Figure 1.Underground hot water storage at the initial state t=0.

Figure 2 .
Figure 2. Left side wall composition used to calculate the current boundary temperature T 1 .

Figure 5 .
Figure 5.Time dependent minimum (orange) and maximum (blue) of the dimensionless stream function ψ.

Figure. 13 .
Figure. 13.Temperature field T COND for solid body cooling.

Table 2 .
Temperature ( o C) and maximum velocity (cm / min) values at different days for Ra = 10 9

Table 3 .
Temperature ( o C) and maximum velocity (cm / min) values at different days for Ra = 10 10