# CFD study of thermocline formation in stratified water storage: Consideration of a second-order Boussinesq approximation to model buoyancy effects and its application to assess the impact of operating conditions

Alexis Ferre, Jérôme Pouvreau, Sylvain Serra, Remi Manceau, Arnaud Bruch

► To cite this version:

Alexis Ferre, Jérôme Pouvreau, Sylvain Serra, Remi Manceau, Arnaud Bruch. CFD study of thermocline formation in stratified water storage: Consideration of a second-order Boussinesq approximation to model buoyancy effects and its application to assess the impact of operating conditions. Proceedings of the 17th International Heat Transfer Conference, IHTC-17 14 – 18 August 2023, Cape Town, South Africa, Aug 2023, Cape Town, South Africa. hal-04145083

HAL Id: hal-04145083

https://hal.science/hal-04145083v1

Submitted on 28 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CFD STUDY OF THERMOCLINE FORMATION IN STRATIFIED WATER STORAGE: CONSIDERATION OF A SECOND-ORDER BOUSSINESQ APPROXIMATION TO MODEL BUOYANCY EFFECTS AND ITS APPLICATION TO ASSESS THE IMPACT OF OPERATING CONDITIONS

Alexis FERRE<fnref n="1" />*, Jérôme POUVREAU<fnref n="1" />, Sylvain SERRA<fnref n="3" /> Rémi MANCEAU<fnref n="2" />, Arnaud BRUCH<fnref n="1" />

## ABSTRACT

Thermal storages are components used in energy systems, such as district heating networks or thermal power plants, in order to decouple the supply of heat from its use. Usage rate of monophasic thermocline storages is highly dependent on the thermal gradient zone inside the fluid, also named *thermocline*. While thermal stratification results of a formation phase followed by a degradation phase, the early stages of thermocline establishment is primarily responsible for its thickness. CFD allow to consider the multiple physical phenomena involved during the thermocline formation, in particular the buoyancy effects. These effects are usually considered by selecting either a variable density with respect to the temperature or a constant one by using the commonly used (first-order) Boussinesq approximation. However, the former approach implies an increased computational cost, and the latter is only valid for an unclear validity range of temperature difference. Hence, this article suggests the use of a second-order Boussinesq approximation, coupled with a RANS turbulence approach, to better account for buoyancy effects in a turbulent water flow submitted to a large temperature differences. CFD results obtained with a quadratic Boussinesq approximation are similar to the one obtained with a variable density but with a computation time reduced by half. This approach is applied to the issue of reducing the thermocline thickness during its creation and the impact of linear flow rate ramps is assessed on both a uniform and initially stratified storage. On an initially cold tank, results show that the longer the ramp time, the thinner the thermocline. In contrast, on the initially stratified tank tested, a gradual injection shows no significant reduction of the thermocline thickness. This can be relevant when performing storage management enhancement.

**KEY WORDS:** Thermal storage, Thermocline, Second-order Boussinesq approximation, Flow rate ramps, Thermal initial state.

## 1. INTRODUCTION

In a monophasic thermocline storage, the usage rate (ratio between the amount of energy effectively stored and the maximum amount of storable energy) is closely dependent on the thermal gradient zone, also named *thermocline*. The more the gradient is spread in the component, the less energy stored can be valorized in the process. In such storages, several physical phenomena are responsible for a thick thermocline: mixing induced by the injection/withdrawal of the fluid at the inlet/outlet, conduction inside the fluid, natural convection during stand-by phases linked to the cooling of the liquid near the wall due to heat losses, and thermal losses among the walls. The first phenomenon is referred as *flow distribution* and is accountable for thermocline formation while the others are responsible for its enlargement. Literature shows that thermocline thickness mainly depends on the formation stage [1].

Various studies have been performed to improve thermocline storages by optimizing its geometry [2] and its tank aspect ratio [3]. However, once it is built, tank structure is fixed and operating conditions

<footnote n="1">¹Univ. Grenoble Alpes, CEA, Liten, DTCH, 38000 Grenoble, France</footnote>

<footnote n="2">²Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, Inria, équipe CAGIRE, LMAP, Pau, France</footnote>

<footnote n="3">³Université de Pau et des Pays de l'Adour, E2S UPPA, LaTEP, Pau, France</footnote>
parameters such as the nominal flow rate and the temperature difference between the inlet and stagnant fluid (Δ*T*) are imposed by the process. Hence, hydraulic distribution devices, which can be subsequently adapted, and the evolution of flow rate before it reaches nominal value, are levers for a better stratification.

The flow rate is linked to the velocity at the inlet/outlet of the tank and is responsible for the turbulent mixing intensity. The greater the mixing, the thicker the thermocline zone. However, if the flow rate is too low, fluid transit time will be substantial and dissipation by conduction might be significant. Sensitivity studies on the flow rate have been performed [4] by varying the nominal value on an initially uniformly cold tank. To the authors' knowledge, no numerical studies on thermocline storages have been carried to assess the impact of the flow rate evolution before reaching its nominal. Moreover, an increase of Δ*T* generates an increase of density difference which stabilizes the flow and ultimately reduces thermocline thickness. But the augmentation of Δ*T* intensifies conduction within the fluid which tends to spread the thermal gradient. Hence, this article intends to study the effects of a gradual flow rate injection with consideration to the initial thermal storage state (uniform or stratified).

CFD is employed because it allows to represent the physical phenomena mentioned above. Because of Δ*T*, buoyancy effects influence the fluid motion. When it comes to water thermocline storages, buoyancy effects are usually considered by choosing either a variable density with respect to temperature [4] or a constant one by using the common first-order Boussinesq approximation (B1) [5]. However, the former approach implies an increased computational cost, and the latter is only valid for an unclear validity range of temperature difference. The B1 approximation assumes a linear density state equation with regards to temperature. To capture the specific density behaviour of some fluids, studies were performed using a non-linear density state equation. For water around 4°C, the objective was to study inversion problem at maximum density temperature with non-linear [6], 2^nd^ order [7] or a 4^th^ order [8] density state equation. In addition, studies have been performed with 2^nd^ order density state equation to model the behaviour of specific fluids such as nanofluids [9] and ferofluids [10] in porous media. To the authors' knowledge, no studies have been performed with a second-order Boussinesq approximation (B2) to capture the behaviour of water with a large temperature difference in a RANS approach. Hence, this article suggests the use of such a method to capture the quadratic behaviour of density on a thermal storage with a Δ*T* of 25°C. A comparison in terms of accuracy and computation time is done between Boussinesq approximations (first and second order) and a variable density. Calculation are based on the case of Zurigat et al. [11] and their experimental data are used for model validation.

## 2. CFD MODEL DETAILS

**2.1 Application Case** Among the data available in literature [4,5,12], the study of Zurigat et al. [11] is chosen because the entire geometry is known (tank, hydraulic distribution system, and sensors), and operating conditions approximate those used by storages in operation. Their case study shows the flow inside a 144.7cm high cylindrical tank with a radius of 20.3cm. Hot water (50.8°C) enters a storage fully filled with cold water (25.9°C), thanks to an impinging jet at the top of the tank (flow is represented on Fig. 1a). The permanent inlet flow rate is 5.92 L/min which corresponds to a vertical fluid velocity inside the storage (*U*<sub>fluid</sub>) of 0.76 mm/s. In this configuration, the Reynolds numbers based on the hydraulic diameter of the inlet pipe and the vessel are 12,816 and 453 respectively. Thus, a coexistence of turbulent and laminar regions is observed and the need to consider turbulence effects is outlined. Instrumentation is positioned at 9 levels (Fig. 1b). At each level, two temperature sensors are located at 5cm and 7cm from the wall.
<figure><img src="image_2.png" /><figcaption>Fig. 1 (a) Sectional view of the central plane of the storage with the flow, the fluid domain and boundary conditions. 2D-axisymetric domain with (b) temperature sensors and (c) generated mesh.</figcaption></figure>

**2.2 Governing Equations** To study the thermocline formation, the following hypotheses are made: flow is assumed to be 2D-axisymmetric; thermal losses, inertia at the walls, transfer by radiation and thermal effects of viscous dissipation are neglected. For a non pressurized water system between 25°C and 51°C, the relative deviation from the average value for molecular viscosity, thermal conductivity and specific heat capacity are respectively -49.3%, 6.00%, 0.01% [13]. Therefore, thermophysical properties are considered variable with the following temperature dependency:

$$ \mu(T) = 2,517 \times 10^{-2} - 1,444 \times 10^{-4}T + 2,101 \times 10^{-7}T^2 \quad (1) $$

$$ \lambda(T) = -0,8492 + 8,046 \times 10^{-3}T - 1,061 \times 10^{5}T^{2} \quad (2) $$

$$ C_p(T) = 1,143 \times 10^{-4} - 65,94T + 0,1991T^{2} - 1,995 \times 10^{-4}T^{3} \quad (3) $$

To simplify the notation, the instantaneous velocity, the average velocity and the fluctuating velocity are noted $u_i^*$, $U_i$ ($U_i = \overline{u_i^*}$) and $u_i$ respectively. The Reynolds decomposition is therefore written $u_i^* = U_i + u_i$. Classical equations of conservation of mass, momentum and energy are solved and expressed below in Cartesian coordinates with Einstein's notation. For turbulence, the k-ω/SST model (Menter, 1994) is retained and equation are described in [14]. To close the equation system, density state equations forms (variable or constant) are discussed in section 3 and the expression of the buoyancy source term $S_b$ in equation (5) is presented in Table 1.

$$ \frac{\partial \rho}{\partial t} + \frac{\partial \rho U_i}{\partial x_i} = 0 \quad (4) $$

$$ \frac{\partial \rho U_i}{\partial t} + \frac{\partial \rho U_i U_j}{\partial x_j} = -\frac{\partial P}{\partial x_i} + \frac{\partial \tau_{ij}}{\partial x_j} + \frac{\partial (-\rho \overline{u_i u_j})}{\partial x_j} + S_b \quad (5) $$

$$ \frac{\partial T}{\partial t} + U_i \frac{\partial T}{\partial x_j} = \frac{1}{\rho} \frac{1}{C_p} \frac{\partial}{\partial x_i} \left( \lambda \frac{\partial T}{\partial x_i} \right) + \frac{\partial}{\partial x_i} \left( \alpha_t \frac{\partial T}{\partial x_i} \right) \quad (6) $$

**2.3 Mesh Generation** Mesh is composed of triangular cells with layers of rectangular cells near the walls with an expansion chosen in order to satisfy the constraint $y^+ \approx 1$ at the wall (Fig. 1c). Rectangular layer thickness is 3mm. The mesh is divided into two zones with a demarcation line chosen at 115cm from the bottom of the tank in order to analyse the thermocline formation while avoiding effects of digital scattering related to mesh size change. By defining a smaller mesh size in the upper part (thermocline formation zone) than in the lower part (thermal front shift), this demarcation allows to refine the zone of interest and to limit the total number of cells in the domain.

With the B1 approximation for modelling buoyancy and constant thermophysical properties, mesh independency study was performed on the mesh sizes $420 \times 10^3$, $680 \times 10^3$, $1200 \times 10^3$ cells. They correspond to a characteristic cell size in the upper part of 0.70mm, 0.50mm and 0.35mm respectively. The size in the lower part remains the same (1.8mm). As the meshes with $680 \times 10^3$ cells and $1200 \times 10^3$ cells did not provide improved results, the mesh with $420 \times 10^3$ cells was adopted for further analysis, maintaining the best trade-off between accuracy and computational resources.
**2.4 Boundary Conditions** The following boundary conditions (Fig. 1a) are applied: inlet fluid temperature is imposed and remains constant ($T_{in} = 50.8^{\circ}C$); the flow regime is considered as fully developed ($I = 5\%$; $D_h = 1,8cm$) and the turbulent variables are calculated according the following equations $k = (3/2)(U_{nom}I)^2$ and $\omega = \sqrt{k}/[C_{\mu}(0,07D_h/C_{\mu}^{3/4}]$; inlet velocity ($U_{in}$) is imposed constant or following a linear ramp profile.

**2.5 Numerical Procedure** The software used for the simulations is Ansys Fluent (2021R2). For the pressure-velocity coupling the pressure-based solver is used with the SIMPLE scheme. A second order upwind type discretization scheme and an implicit first-order temporal discretization are used to solve the transport equations. In order to reduce computational time, the adaptive error-based time step is used. The initial time step is equal to $1 \times 10^{-3}$s. To limit and control numerical error, the maximum time step is limited such as not to exceed a maximum CFL number of 100 in the domain. The maximum time step is equal to $1 \times 10^{-2}$ s.

## 3. ACCOUNTING FOR BUOYANCY EFFECTS

**3.1 Variable Density** To model buoyancy in CFD, two approaches can be thought out: considering a temperature dependent density or a constant one. Between 25°C and 51°C, non-pressurized water density variations are accurately represented by a quadratic behaviour. With the Boussinesq approximation, density can be expressed with regard to a reference temperature $T_0$ which reveals the thermal expansion coefficient $\beta_1$ obtained at $T_0$:

$$ \rho(T) = \rho_0 - \beta_1\rho_0(T - T_0) - \beta_2\rho_0(T - T_0)^2 \quad (7) $$

Using a temperature variable density allows to obtain buoyancy term source $S_b = [\rho(T) - \rho_0]g_i$ without assumptions. However, this method requires more computation resources.

### 3.2 Constant Density

**1st order Boussinesq approximation.** The second way to model buoyancy is by considering a constant density and using the Oberbeck-Boussinesq approximation [15,16]. This approximation neglects variation of density except when it is directly responsible for buoyancy. In addition, the common approximation relies on the following assumptions: small temperature differences, a linear density state equation, negligible viscous heat dissipation, constant thermophysical properties and small hydrostatic pressure variations. In the common Boussinesq approximation, a Taylor series expansion of density is performed at the reference temperature (8) and only the first order term is kept.

$$ \rho(T) = \rho_0 + \rho'(T_0)(T - T_0) + \frac{1}{2}\rho''(T_0)(T - T_0)^2 \quad (8) $$

While this method reduces the computational cost, it assumes a linear buoyant source term which ignores the quadratic correction. This is a reason why the B1 approximation is only valid for a small temperature difference. Gray and Giorgini [17] studied the temperature validity range and concluded that the temperature difference for problems involving air and water are less than 28.6°C and 1.25°C respectively. Ferziger and Perić [18] stated that the use of the Boussinesq approximation ‘introduces errors of the order of 1% if the temperature differences are below e.g. 2°C for water and 15°C for air’. Thus, studies tried to improve this approach.

**2nd order Boussinesq approximation.** In their 2021 review, Mayeli and Sheard [19], suggest a classification of the different approaches of the Boussinesq approximation from the original onto modified versions which aim to improve it. In the present study, when the temperature difference becomes large enough, density and other thermophysical properties are no longer linear with temperature and higher term cannot be neglected in the density state equation. Hence, the second order
term of the Taylor series expansion should be considered to capture the quadratic density behaviour. Density state equation and buoyancy terms are synthesized in **Table 1**.

<table><caption>**Table 1** Expressions of the studied buoyancy models</caption><thead><tr><th>Model</th><th>Density</th><th>Buoyancy term<br>S<sub>b</sub> = -[ρ(T) - ρ<sub>0</sub>]g</th></tr></thead><tbody><tr><td>ρ(T)</td><td>ρ<sub>0</sub> - β<sub>1</sub>ρ<sub>0</sub>(T - T<sub>0</sub>) - β<sub>2</sub>ρ<sub>0</sub>(T - T<sub>0</sub>)<sup>2</sup></td><td>β<sub>1</sub>ρ<sub>0</sub>(T - T<sub>0</sub>)g + β<sub>2</sub>ρ<sub>0</sub>(T - T<sub>0</sub>)<sup>2</sup>g</td></tr><tr><td>B1</td><td>ρ<sub>0</sub></td><td>β<sub>1</sub>ρ<sub>0</sub>(T - T<sub>0</sub>)g</td></tr><tr><td>B2</td><td>ρ<sub>0</sub></td><td>β<sub>1</sub>ρ<sub>0</sub>(T - T<sub>0</sub>)g + β<sub>2</sub>ρ<sub>0</sub>(T - T<sub>0</sub>)<sup>2</sup>g</td></tr></tbody></table>

### 3.3 Buoyancy Models Comparison

Temperature fields obtained with the three buoyancy models mentioned above are compared at 60s to focus on the thermocline formation (Fig. 2). Temperature fields with a variable density is used as a reference to plot the difference between the B1 and B2 approximations. The ones obtained with a varying density and the B2 approximation are almost identical (min = -1.08°C, max = 2.19°C). In contrast, the temperature fields with B1 approximation show significant differences (min = -9.06°C; max = 6.03°C). The velocity magnitude, turbulent kinetic energy, specific dissipation rate have been compared and show comparable deviation as the one obtained with temperature (not shown here).

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Comparison of temperature fields for ρ(T), B2 and B1 buoyancy approaches (t=60s)</figcaption></figure>

The computational time between a variable density, the B2 and B1 approximations are compared. The time obtained is based on 60s of physical simulation. The computational time obtained with a variable density, the B1 and B2 approximations are 9792s, 4236s and 4550s respectively. The use of the B2 approximation allows a computation time reduction of 54% compared to a variable density. Since the B2 approximation provides an accurate representation of buoyancy effects for a large temperature range with a significantly reduced computational time, it is used in the rest of the present study.
## 4. COMPARISON WITH EXPERIMENTS

In order to validate the CFD model, a comparison between experimental and numerical results is performed in this section. To reproduce the experimental settings, the numerical setup meets the following requirements: the calculated temperature is the average temperature between $TC_{5cm}$ and $TC_{7cm}$; the inlet flow instantaneously reaches its nominal value. Experimentally, the tank is supposed to be initially filled with cold water. However, experimental data (Fig. 3) show that the initial temperature at $TC_{9}$ is higher than the initially announced temperature (27.1 °C instead of 25.9°C) and the temperature immediately starts rising when fluid is injected. Thus, the storage may not be initially uniformly cold. Hence, an initially stratified storage is considered. To thermally initialize the storage, a temperature gradient is rebuilt from the experimental values. The obtained temperature field function is $T = T_{cold} + [(z - c)/a]^{1/b}$ with $a = 1.15 \ m.K^{-b}$, $b = 0.05$ and $c = 0.12 \ m$. The time evolution of the temperature is in Fig. 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Comparison between numerical and experimental results</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 Thermal initial states: (a) uniform and (b) stratified</figcaption></figure>

With a uniform initial state (Fig. 4a), the numerical temperature evolution at level 9 exhibits irregularities that cannot be seen in the experimental data (Fig. 3). These irregularities suggest the existence of eddies in the upper part of the tank. Moreover, even if the configuration is able to capture the shape of the thermal gradient, the rise of the numerically-obtained temperature is slower than the experimental one. With an initial gradient (Fig. 4b), experimental and numerical results show a better agreement with a decrease of the temperature rise delay. The existence of a gradient must be accounted for and this approach is applied to study the impact of the operating conditions on the thermocline formation in the next section.

## 5. IMPACT OF A GRADUAL FLOW RATE INJECTION

### 5.1. Thermocline Thickness
To quantify operating condition effects, the thermocline thickness based on temperature limits [20] is selected. This criterion relies on the admissible hot and cold limits of temperature and focuses on thermocline edges (Fig. 6). Usually, the industrial process in which the storage is included defines the fraction of the tolerable temperature difference. Temperature limits are assessed with the dimensionless temperature $\theta$.

$$ \theta(z) = \frac{T(z) - T_{\text{cold}}}{T_{\text{hot}} - T_{\text{cold}}} \quad (9) $$

$$ \delta = U_{\text{fluid}}\Delta t \quad (10) $$

Denominator represents the maximal temperature difference of the storage. In this study, the chosen dimensionless temperature range is $0.15 \le \theta \le 0.85$. In other words, thermocline represents 70% of the maximal temperature difference symmetrically distributed. With a dimensional temperature, the thermocline thickness corresponds to $29.6^{\circ}\text{C} \le T(z) \le 47.1^{\circ}\text{C}$. Experimentally, the temperature
evolution at a given point is usually obtained and an ideal plug flow is assumed to reconstruct thermocline the thickness (10). Here, the temperature evolution at $TC7_{7cm}$ is used.

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Thermocline thickness under the dimensionless temperature</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Impact of flow rate ramps and initial conditions on thermocline thickness</figcaption></figure>

**5.2. Initially Uniform Storage** Since it is the usual configuration of experimental characterization tests, the impact of a gradual injection on a homogenously initially cold tank is assessed in this section. Gradual injection is done according to a linear flow rate ramp. It means that flow rate needs $t_{ramp}$ seconds to go from zero to the nominal value. Five values of $t_{ramps}$ are tested: 0.1s, 30s, 90s, 120s and 240s. Because the flow rate is linearly increasing before reaching its permanent value, variables are compared at same amount of injected volume to compare at the same amount of stored energy. Qualitatively, temperature evolution shows that the thermocline thickness significantly decreases with the increase of $t_{ramps}$ and an asymptote seems to be reached at 240s (Fig. 7a). This is quantitatively confirmed because the use of a 240s linear ramp reduces the thickness by 73% (Fig. 6). A ramp time of 240s represents 12.6% of the total filling time. To understand the thermocline creation process, temperature fields and velocity vectors are plotted for different values of injected volume for two ramp times: 0.1s (abrupt injection) and 120s (gradual injection) (Fig. 8).

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Temperature evolution for TCs located at 7cm from the storage wall with an initially (a) homogenously cold and (b) stratified storage</figcaption></figure>

With an abrupt injection, the incoming fluid immediately impacts the side wall of the storage and creates a large eddy ($V_1$). Levers and Lin [21] also show in their CFD study that after hitting the reservoir wall, the jets will spread over a large area and significantly destroy the level of thermal stratification. And, the more the flow rate increases, the more the stratification is affected. In the studied case, as incoming fluid gets inside the storage, further mixing is induced ($V_2$) and the thermocline thickens. Velocity vectors show that $\Delta T$ affects eddy's shape and prevents its formation. At $V_4$, the thermocline is almost formed. At $V_5$, it reaches its final shape.
With a gradual injection, incoming water allows to immediately create a thin thermocline ($V_1$). The phenomena governing the generation of vortices could explain this result: the vortices is only generated by the setting in motion of the initial fluid than by the turbulence phenomena [22]. So a gradual injection generate smaller vortices than an injection with a flow rate step, and therefore limit mixing. Another perspective is to explain gradient formation with the density (or gravity) current theory. This current is an horizontal flow in a gravitational field that is driven by a density difference in a fluid and is constrained to flow horizontally (e.g. hot air on a ceiling). In 1986, Yoo et al. [1] studied thermocline formation and concluded that good stratification will appear when fluid exiting the hydraulic distributor forms a gravity current. It has been observed that the first eddy triggers a maximum penetration height ($h_{max}$) when impacting the wall (Fig. 8). This height is mainly responsible for the density current's thickness [23]. Applied here, the more gradual the injection, the thinner the penetration height and the less intense the mixing. The creation of a thin thermocline during the ramp enables to limit the mixing by strongly preventing the eddy to be fully formed ($V_2$). Afterwards, thermocline is formed ($V_3$) and thickens due to thermal conduction ($V_4$). However, during operating cycles of an industrial process, thermocline storages rarely discharge the entirety of the thermocline. Hence, charging phases are tested on an initially stratified storage.

<figure><img src="image_9.png" /><figcaption>Fig. 8 Temperature field and velocity vector during thermocline formation for a 0.1s and 120s flow rate injection</figcaption></figure>

**5.3. Initially Stratified Storage** Impact of initial thermal stratification is explored in this section. Qualitatively, the use of ramps does not have a significant impact on temperature evolution (Fig. 7). Quantitatively, the thermocline thickness remains constant regardless of the ramp time (Fig. 6). It can be relevant for storage operating conditions management to perform at high flow rate when the storage is already stratified. It allows to deliver a high nominal power. Therefore, an improve management storage can be to primary create a thin thermocline with a gradual injection on an initially cold tank. Then, once the tank is stratified, the nominal flow rate can be employed in the industrial process.
# 6. CONCLUSION

A CFD study on thermocline formation have been performed. Three different ways of considering buoyancy effects have been compared. Results show that the use of Boussinesq approximation with a quadratic density state equation allow to correctly represent buoyancy effects while reducing the computational time by a factor 2 compared to a variable density calculation. In addition, effect of linear flow rate ramps have been numerically tested on both an initially uniformly cold and stratified tank. If the storage is fully filled with cold water, ramps allow to create a thin thermocline because it smoothly generates a strong thermal gradient which prevent eddies to fully develop. The formation of a thin thermocline during the early stages can be explained by density current theory. With a stratified storage, ramps have no significant impact. To further study the effect of initial stratification, another initial state will be tested in the future. Present results suggest that with a proper storage management, the thermocline thickness can be reduced. Finally, while the effects of ramps and initial thermal states have been numerically investigated, an experimental verification is necessary. Hence, an experimental apparatus is being dimensioned within the scope of the project.

# NOMENCLATURE

## Latin symbols

$C_p$ specific heat capacity [J.K⁻¹.kg⁻¹]

$C_μ$ dimensionless coefficient, [-]

$g$ gravity acceleration [m.s⁻²]

$I$ turbulent intensity [-]

$k$ turbulent kinetic energy [J.kg⁻¹]

$P$ pressure [Pa]

$S_{ij}$ mean strain rate tensor, [s⁻¹]

$t$ time [s]

$T$ temperature [K]

$U$ mean velocity [m.s⁻¹]

$u$ fluctuating velocity [m.s⁻¹]

$u^*$ instantaneous velocity [m.s⁻¹]

$\overline{u_i u_j}$ Reynolds stress tensor, [m².s⁻²]

$y^+$ dimensionless wall distance [-]

$z$ axial coordinate [m]

## Greek symbols

$α$ diffusivity [m².s⁻¹]

$β$ thermal expansion coefficient, [-]

$δ_{ij}$ Kronecker symbol [-]

$λ$ thermal conductivity [W.m⁻¹.K⁻¹]

$μ$ dynamic viscosity [Pa.s]

$ν$ kinematic viscosity [m².s⁻¹]

$ρ$ density [kg.m⁻³]

$τ_{ij}$ viscous shear stress, [Pa]

$τ_{ij,t}$ turbulent shear stress, [Pa]

$ω$ specific dissipation rate [s⁻¹]

## Indices and exponents

$0$ reference

$i,j$ cartesian directions

$rel$ relative

$t$ turbulent

## Abbreviations

$B1$ 1st order Boussinesq

$B2$ 2nd order Boussinesq

$CFD$ Computational Fluid Dynamics

$IC$ Initial Conditions

$RANS$ Reynolds-Averaged Navier-Stokes

$Δ$ Variation

# REFERENCES

[1] J. Yoo, M. W. Wildin, and C. R. Truman, 'Initial formation of a thermocline in stratified thermal storage tanks', ASHRAE Trans U. S., vol. 92:2A, Art. no. CONF-8606125-, Jan. 1986, Accessed: May 19, 2021. [Online]. Available: https://www.osti.gov/biblio/6960784-initial-formation-thermocline-stratified-thermal-storage-tanks

[2] Z. Yang, H. Chen, L. Wang, Y. Sheng, and Y. Wang, 'Comparative study of the influences of different water tank shapes on thermal energy storage capacity and thermal stratification', Renew. Energy, vol. 85, pp. 31–44, Jan. 2016, doi: 10.1016/j.renene.2015.06.016.

[3] J. E. B. Nelson, A. R. Balakrishnan, and S. S. Murthy, 'Experiments on stratified chilled-water tanks', Int. J. Refrig., p. 19, 1999.
[4] Y. P. Chandra and T. Matuska, 'Numerical prediction of the stratification performance in domestic hot water storage tanks', *Renew. Energy*, vol. 154, pp. 1165-1179, Jul. 2020, doi: 10.1016/j.renene.2020.03.090.

[5] L. J. Shah and S. Furbo, 'Entrance effects in solar storage tanks', *Sol. Energy*, vol. 75, no. 4, Art. no. 4, Oct. 2003, doi: 10.1016/j.solener.2003.04.002.

[6] B. Gebhart and J. C. Mollendorf, 'A new density relation for pure and saline water', *Deep Sea Res.*, vol. 24, no. 9, pp. 831-848, Sep. 1977, doi: 10.1016/0146-6291(77)90475-1.

[7] C. V. Raghavarao and Y. V. S. S. Sanyasiraju, 'Natural convection heat transfer of cold water between concentric cylinders for high rayleigh numbers—a numerical study', *Int. J. Eng. Sci.*, vol. 32, no. 9, pp. 1437-1450, Sep. 1994, doi: 10.1016/0020-7225(94)90122-8.

[8] P. Vasseur, L. Robillard, and B. Chandra Shekar, 'Natural Convection Heat Transfer of Water Within a Horizontal Cylindrical Annulus With Density Inversion Effects', *J. Heat Transf.*, vol. 105, no. 1, pp. 117-123, Feb. 1983, doi: 10.1115/1.3245529.

[9] K. Thriveni and B. Mahanthesh, 'Significance of variable fluid properties on hybrid nanoliquid flow in a micro-annulus with quadratic convection and quadratic thermal radiation: Response surface methodology', *Int. Commun. Heat Mass Transf.*, vol. 124, p. 105264, May 2021, doi: 10.1016/j.icheatmasstransfer.2021.105264.

[10] H. M. Elshehabey, Z. Raizah, H. F. Öztop, and S. E. Ahmed, 'MHD natural convective flow of Fe3O4-H2O ferrofluids in an inclined partial open complex-wavy-walls ringed enclosures using non-linear Boussinesq approximation', *Int. J. Mech. Sci.*, vol. 170, p. 105352, Mar. 2020, doi: 10.1016/j.ijmecsci.2019.105352.

[11] Y. H. Zurigat, P. R. Liche, and A. J. Ghajar, 'Influence of inlet geometry on mixing in thermocline thermal energy storage', *Int. J. Heat Mass Transf.*, vol. 34, no. 1, pp. 115-125, Jan. 1991, doi: 10.1016/0017-9310(91)90179-I.

[12] A. Zachár, I. Farkas, and F. Szlivka, 'Numerical analyses of the impact of plates for thermal stratification inside a storage tank with upper and lower inlet flows', *Sol. Energy*, vol. 74, no. 4, Art. no. 4, Apr. 2003, doi: 10.1016/S0038-092X(03)00188-9.

[13] 'National Institute of Standards and Technology', NIST. https://www.nist.gov/ (accessed Nov. 23, 2022).

[14] F. R. Menter, 'Two-equation eddy-viscosity turbulence models for engineering applications', *AIAA J.*, vol. 32, no. 8, pp. 1598-1605, Aug. 1994, doi: 10.2514/3.12149.

[15] A. Oberbeck, 'Ueber die Wärmeleitung der Flüssigkeiten bei Berücksichtigung der Strömungen infolge von Temperaturdifferenzen', *Ann. Phys.*, vol. 243, no. 6, pp. 271-292, 1879, doi: 10.1002/andp.18792430606.

[16] J. (1842-1929) A. du texte Boussinesq, *Théorie analytique de la chaleur mise en harmonie avec la thermodynamique et avec la théorie mécanique de la lumière*, Gauthier-Villars., vol. 2, 2 vols. Paris, 1903. Accessed: Dec. 03, 2022. [Online]. Available: https://gallica.bnf.fr/ark:/12148/bpt6k61635r

[17] D. D. Gray and A. Giorgini, 'The validity of the boussinesq approximation for liquids and gases', *Int. J. Heat Mass Transf.*, vol. 19, no. 5, Art. no. 5, May 1976, doi: 10.1016/0017-9310(76)90168-X.

[18] J. H. Ferziger and M. Perić, *Computational Methods for Fluid Dynamics*. Berlin, Heidelberg: Springer Berlin Heidelberg, 2002. doi: 10.1007/978-3-642-56026-2.

[19] P. Mayeli and G. J. Sheard, 'Buoyancy-driven flows beyond the Boussinesq approximation: A brief review', *Int. Commun. Heat Mass Transf.*, vol. 125, p. 105316, Jun. 2021, doi: 10.1016/j.icheatmasstransfer.2021.105316.

[20] A. Musser and W. P. Bahnfleth, 'Parametric Study of Charging Inlet Diffuser Performance in Stratified Chilled Water Storage Tanks with Radial Diffusers: Part 1-Model Development and Validation', *HVACR Res.*, vol. 7, no. 1, pp. 31-49, Jan. 2001, doi: 10.1080/10789669.2001.10391428.

[21] S. Ievers and W. Lin, 'Numerical simulation of three-dimensional flow dynamics in a hot water storage tank', *Appl. Energy*, vol. 86, no. 12, Art. no. 12, Dec. 2009, doi: 10.1016/j.apenergy.2009.04.010.

[22] S. Carpy and R. Manceau, 'Turbulence modelling of statistically periodic flows: Synthetic jet into quiescent air', *Int. J. Heat Fluid Flow*, vol. 27, no. 5, pp. 756-767, Oct. 2006, doi: 10.1016/j.ijheatfluidflow.2006.04.002.

[23] E. Kaloudis, D. G. E. Grigoriadis, and E. Papanicolaou, 'Numerical simulations of constant-influx gravity currents in confined spaces: Application to thermal storage tanks', *Int. J. Therm. Sci.*, vol. 108, pp. 1-16, Oct. 2016, doi: 10.1016/j.ijthermalsci.2016.04.018.