Numerical Simulation Study of Thermal Performance in Hot Water Storage Tanks with External and Internal Heat Exchangers

Abstract

This paper presents a numerical analysis of two hot water storage tank configurations—one equipped with an external heat exchanger (Tank-1) and the other with an internal heat exchanger (Tank-2). The objective is to evaluate and compare their thermal performance during charging and discharging processes. The numerical model is developed by solving a system of ordinary differential equations using the 4th-order Runge–Kutta method, implemented in the Python programming language. The results indicate that Tank-1 demonstrated a higher charging efficiency of 94.6%, achieving full charge in approximately 2 h and 20 min. In comparison, Tank-2 required 3 h and 47 min to reach full charge, with a charging efficiency of 85.9%. During discharge, both configurations exhibited similar behavior, with an efficiency of 13.63% over approximately 33 min. The analysis showed that the external heat exchanger configuration led to more effective thermal stratification, supported by the Richardson number analysis, which indicated a significant effect of buoyancy during charging. This design advantage makes Tank-1 particularly suitable for applications requiring rapid heating and minimal heat loss, such as in cold climates or intermittent demand systems. The numerical model demonstrated reliable predictive accuracy, achieving an RMSE of 6.1% for the charging process and 6.8% for the discharging process, thereby validating the model’s reliability. These findings highlight the superior performance of the external heat exchanger configuration for fast and efficient energy storage, particularly for applications in cold climates.

1. Introduction

Thermal energy storage (TES) systems are essential for improving energy efficiency and enabling the integration of renewable energy sources, especially under variable climate conditions and fluctuating energy demands. These systems enable efficient storage of excess energy for later use, reducing dependency on fossil fuels and improving grid stability.

TES is crucial for balancing the intermittent nature of renewable energy sources [1]. A widely implemented TES method is sensible heat storage (SHS), which stores thermal energy by raising the temperature of a storage medium [2]. Water storage tanks (WSTs) are commonly used for SHS due to water’s high specific heat capacity and cost-effectiveness. However, the energy density of SHS is relatively low, requiring larger storage volumes to store substantial energy. To address this limitation, latent heat storage (LHS) materials, such as phase change materials (PCMs) [3], can be integrated into SHS systems. This combination of SHS with LHS, using materials like PCMs, offers substantial advantages by enhancing both thermal stratification [4] and energy density, making it suitable for diverse practical applications, particularly in residential heating, solar thermal systems, and industrial energy storage [5].

TES tanks are generally more cost-effective and easier to implement than solid-state heat accumulators, particularly due to the widespread use of water, which has a high specific heat capacity and is readily available [1,2]. These WSTs are flexible in capacity, scalability, and maintenance, making them suitable for various applications. In contrast, solid-state heat accumulators, which use specialized materials like molten salts or metals, offer higher energy density but come with higher costs and more complex handling requirements, making them less practical for many general energy storage needs [6].

Depending on the heat source and climatic conditions, a hot WST can be connected to the heat source in various ways. For instance, in hot climate regions, solar water heaters are typically of the thermosiphon type [7], where a solar thermal collector is connected to the tank without the need for a circulation pump or intermediate heat exchangers. In this configuration, water is used as the heat transfer fluid (HTF) in both the tank and the solar collector. In contrast, in cold climate regions, an antifreeze fluid is used in the solar collector, which circulates with the help of a pump [8]. This HTF then heats the water in the tank using an intermediate heat exchanger. This heat exchanger can be submerged inside the tank and installed externally to prevent mixing between the antifreeze fluid and the water. The present study focuses on comparing these two connection configurations to numerically assess the efficiency of TES. Such tanks can be utilized in conjunction with solar collectors, heat pumps, and LHS accumulators filled with PCMs. The objective is to evaluate and compare the thermal performance of two hot water storage tanks, one with an internal heat exchanger and the other with an external heat exchanger configuration, during the charging and discharging processes.

Osiadacz et al. [9] have presented a comparison of three one-dimensional models. In the first case, the convection and conduction equation (Fourier–Kirchhoff equation) is used. An implicit discretization scheme was employed to solve the system of equations, and the solution was obtained using a matrix method. The second method involved using the standard convection equation, with a forward difference combined with a central difference discretization. The resulting scheme is explicit and does not require an additional solver for the solution. In the third case, a one-dimensional thermal balance equation was considered, using a forward difference discretization scheme, resulting in an explicit Euler equation. The comparison of all three mathematical models based on water layer temperatures showed only minor differences between them. Oliveski et al. [10] compared two methods for modeling a hot WST: a two-dimensional model in cylindrical coordinates and three one-dimensional thermal balance models combined with artificial intelligence (AI). The two-dimensional model utilized conservation equations for mass, momentum, and energy, employing the finite volume method along with the Launder and Sharma turbulence model. The one-dimensional models included a multi-node inversion method, an averaged inversion with AI, and a standard multi-layer approach. The two-dimensional model showed good agreement with experiments, while the one-dimensional models, especially those with inversion, provided comparable accuracy with significantly faster computation. De la Cruz-Loredo et al. [11] employed a one-dimensional energy conservation model, solved using hybrid continuous-discrete systems. The calculations were performed in MATLAB using its built-in methods. Alongside the modeling, experimental validation was also conducted. Simulations were carried out with divisions of 12 and 60 layers, and it was found that significant errors occurred in the 60-layer simulation with this one-dimensional model. The best agreement with experimental data was achieved with the 12-layer division, with an RMSE not exceeding 7% in both cases. Mongibello et al. [12] presented a one-dimensional modeling of the system along with experimental validation. The study involved a tank with an IHX, and the explicit Euler numerical method was used to solve the one-dimensional thermal balance equations, incorporating various accepted correlations. During validation, the RMSE did not exceed 2% for each type of correlation. Cruickshank and Harrison [13] examined a tank with an external heat exchanger and direct water feed into the tank. The modeling was based on one-dimensional energy balance equations, using a 9-layer division of the fluid in the tank. The study focused solely on the analysis of thermal losses. The one-dimensional model showed good agreement with experimental data regarding heat losses. Mongibello et al. [14] presented a model of a solar water heater with a tank containing an IHX, simulated using one-dimensional energy balance equations. The calculation algorithm for the water in the tank proceeds from bottom to top to account for the upward flow of hot water. The code was implemented in MATLAB, and no comparisons with experimental data were conducted.

Aguilar et al. [15] presented a model for a hot WST alongside a heat pump. From the conservation of energy and momentum equations, they derived a one-dimensional model discretized using the finite volume method. Various correlations were also added to describe the mixing of fluid streams. The model showed good accuracy in predicting stratification in the tank, with a deviation of 2.6 °C when compared to experimental data. Kleinbach et al. [16] used TRNSYS 13.1 software for modeling in their study and described the mathematical model implemented in the program. Experimental validation was conducted, with an error not exceeding 13%. The study examined the optimal number of layers and considered two types of modeling approaches: a multi-node model and a plug flow model. The maximum number of segments for the plug flow model was 50. The primary difference between the two models is that in the plug flow model, segments can change in size or increase/decrease in number (dynamic grid). The concept presented in the paper is that the fluid can move both upwards and downwards if the incoming fluid temperature is lower than the temperature of the upper layer. When comparing the two approaches, the plug flow model shows higher efficiency. However, both models tend to overestimate the temperature. The use of the plug flow model is justified if there are variable inlet positions. A disadvantage of the plug flow model is its high sensitivity to the time step. Rahman et al. [17] presented a one-dimensional model of a hot WST, with 10 fluid layers, based on energy balance equations. Correlation coefficients were also used to describe turbulent and laminar flow conditions, along with considering the effects of friction. The equations are nearly identical to those in our model, with the main differences being in the correlation coefficients and specific equations used. The one-dimensional model was compared with simulations conducted in COMSOL 5.0, and the RMSE relative to the COMSOL simulation was 1.01 K for an 1800-s simulation. Raccanello et al. [18] analyzed the modeling of direct and indirect heating of a hot water tank, with and without PCM materials. One-dimensional, quasi-one-dimensional, and zero-dimensional models were used. The zero-dimensional model represented a fully mixed tank, the quasi-one-dimensional model was used for moving boundaries with two layers, and the one-dimensional model was a multi-node model with fixed boundaries, all based on energy balance equations. Correlation coefficients applicable across multiple methods were discussed. The one-dimensional models were experimentally validated for both charging and discharging, with an error of less than 9%. A comparison of the three methods showed that the zero-dimensional and quasi-one-dimensional models had similar limitations, whereas the one-dimensional model proved to be the most accurate. Rahman et al. [19] presented a model of a hot WST with a dual heat exchanger setup—one circulating hot water and the other cold water. One-dimensional energy balance equations were used for the modeling, with equations like those described earlier [17] and to our model. Numerical calculations were validated against experimental data, with the tank stratified into 10 layers. The RMSE did not exceed 2%. Chandra and Matuska [4] presented a review on thermal stratification in domestic hot water storage tanks. The study employed three modeling methods: the multi-node method, the jet entrainment model, and the plug flow method. All methods are one-dimensional and based on energy balance equations, focusing on the movement of the thermocline. No experimental comparisons were included. The paper also examined various tank configurations and their optimization for maintaining stratification, outlining the advantages and disadvantages of the models, as well as explaining the causes of errors within them. As a review article, no specific numerical results were provided. Mongibello et al. [20], in another study, presented calculations for a solar water heater system using TRNSYS 17 and MATLAB. A mathematical model implemented in the software was provided. The MATLAB code detailed the one-dimensional energy balance equations for each component, solved using the Euler numerical method. A comparison between the MATLAB code and TRNSYS showed a maximum difference of 10%. The tank featured an IHX circulating a 30% propylene glycol-water solution. The two models showed similarities, though differences were observed due to the different computational algorithms and mathematical models used.

These studies highlight the strengths and limitations of various models, particularly in achieving accurate thermal stratification and computational efficiency. However, there remains a gap in understanding the comparative efficiency of TES tank configurations (external vs. internal heat exchangers) under varied environmental conditions, particularly with advanced modeling techniques.

This study presents a distinctive approach to one-dimensional TES modeling by employing the Runge–Kutta method along with a fixed-position adaptive grid—an innovation not seen in the prior literature. Additionally, it incorporates specialized correlation equations for Nusselt number calculation, which further enhances the model’s accuracy. The fluid is divided into 20 layers, enabling detailed stratification analysis without accounting for downward fluid movement. This model balances computational efficiency and precision, producing RMSE and deviation results comparable to those of prior studies. A Richardson number analysis was conducted, and thermal efficiency was assessed for both tank configurations, each with identical physical and geometric properties.

2. Configurations of Hot Water Storage Tanks

This study examines two configurations for connecting heat exchangers in a hot WST. A brazed plate heat exchanger (PHX) is employed as an intermediate heating device between the heat source and the storage tank in the first configuration. The intermediate heat transfer fluid (HTF), typically an antifreeze solution, circulates between the heat source—such as solar thermal collectors or heat pumps—and the PHX. The water within the tank circulates between the brazed heat exchanger and the tank itself. Notably, this configuration does not include a heat exchanger within the tank. Heated water from the PHX is introduced at the top region of the tank, while the cooler return water to the heat exchanger is drawn from the bottom region of the tank. Figure 1a illustrates the first configuration, where the heat exchanger is positioned outside the tank. This configuration is referred to as Tank-1 in the following discussion.

Tank-1 is advantageous in applications where ease of maintenance and flexibility in heat exchanger positioning are priorities. This external configuration allows for simpler repairs or modifications without disrupting the tank’s internal structure, making it suitable for retrofitting in existing systems.

For comparison, the second configuration, shown in Figure 1b, incorporates a submerged heat exchanger. In this setup, a non-freezing fluid circulates within the immersed heat exchanger (IHX), serving as the HTF between the heat source—such as a solar collector or heat pump—and the hot WST. The heat exchanger is typically positioned at the bottom of the tank to promote even water heating and minimize thermal stratification. This configuration is designated as Tank-2. In both configurations, cold mains water enters the tank from the bottom, while the hot water outlet for use is located at the top of the tank.

Tank-2’s internal heat exchanger setup promotes higher thermal efficiency by minimizing thermal losses often associated with external heat transfer. This setup is particularly suitable for applications in space-constrained environments, where compactness and efficient heat retention are essential. By positioning the IHX at the tank’s bottom, this configuration optimizes natural convection, supporting a more stable temperature gradient that is beneficial in systems requiring sustained thermal stratification.

To facilitate a fair comparison between the two configurations, both tanks were modeled with identical volumes and geometric properties. The operating conditions were also standardized to isolate the performance impact of each heat exchanger setup. The following fixed parameters were selected for the charging process: hot water at a constant temperature of 45 °C is supplied at a steady flow rate of 1.67 × 10⁻² m³/s for both tanks, with and without the IHX. During the simulation in Tank-2, water was used as the HTF in the IHX. The discharge mode is assumed to be identical for both tanks, so no differences are expected between the two types. The volume occupied by the heat exchanger in the second tank is considered negligible, with the tank containing 151 L of water. The volume of the IHX is 0.0075 m³ or 7.5 L. This is less than 5% of the total volume, resulting in minimal differences in the numerical calculation outcomes for 151 L and 143.5 L. The charging criterion is defined as reaching an average temperature of 40.5 °C, corresponding to 90% of the supplied hot water temperature. The discharge criterion is achieving an average water temperature of 16.5 °C, representing 90% of the incoming cold-water temperature of 15 °C. The tank dimensions are 1.3 m in height and 0.4 m in diameter.

While this comparative analysis provides insights into the performance of both configurations under controlled conditions, it does not account for variables such as fluctuating heat source temperatures or variations in flow rate. In both configurations, circulation pump effects are excluded to focus solely on the natural convection dynamics within each tank. This assumption simplifies the model and highlights the inherent stratification properties of each setup, independent of forced circulation influences.

3. Heat Balance Model Equations

For Tank-1, where the water supply is direct, the energy balance equation is expressed as follows:

$$C_w{m}_{w_i}{\displaystyle \frac{d{T}_{w_i}}{dt}}={\dot{Q}}_{c,w_i}+{\dot{Q}}_{cf, w_i}-{\dot{Q}}_{{loss}_i}-{\dot{Q}}_{hwc},$$

(1)

where ${\dot{Q}}_{c,w_i}$ is the convective heat transfer rate between the liquid layers$, {\dot{Q}}_{cf,w_i}$ is the forced convective heat transfer rate due to direct addition and withdrawal of water, ${\dot{Q}}_{{loss}_i}$ is the heat loss rate to ambient, ${\dot{Q}}_{hwc}$ is the forced convective heat transfer rate due to hot water intake and cold water addition. The index $i$ denotes the layers of liquid in the tank and take values from 1 to 20, respectively.

The convective heat transfer rate between layers is calculated as:

$${\dot{Q}}_{c,w_i}=\left({\displaystyle \frac{\left(T_{w_{i+1}}-T_{w_i}\right)k_w{\mathsf{\Delta }}_{i+1}}{L_w}}+{\displaystyle \frac{\left(T_{w_{i-1}}-T_{w_i}\right)k_w{\mathsf{\Delta }}_{i-1}}{L_w}}\right)A_{w,}$$

(2)

where $A_w$ is the cross-sectional area of the tank, $\mathsf{\Delta }$ is the Adaptive Grid multiplier, $L_w$ is the thickness of the liquid layers, $k_w$ is the thermal conductivity coefficient of water.

The heat transfer rate ${\dot{Q}}_{c,w_i}$ between layers utilizes an Adaptive Grid approach presented by Powell and Edgar [21], which enhances heat transfer by increasing the mixing between layers whenever a temperature difference exists. The $\mathsf{\Delta }$ multiplier, calibrated at a value of 400 based on experimental results, enables more efficient stratification in the tank by accommodating natural convection effects more realistically.

The term for forced convective heat transfer rate is calculated using by:

$${\dot{Q}}_{cf,w_i}=\dot{m}{C}_w\left(T_{w_{i-1}}-T_{w_i}\right),$$

(3)

where $\dot{m}$ is the flow rate of incoming water.

The term ${\dot{Q}}_{cf,w_i}$ captures the rate of temperature change in the upper layers due to incoming hot water (45 °C), promoting stratification and supporting thermal efficiency by counteracting heat losses to the ambient.

The heat loss rate is defined as:

$${\dot{Q}}_{{loss}_i}={\displaystyle \frac{A_{{wall}_i}{k}_{ti}\left(T_{w_i}-293.15\right)}{L_{ti}}},$$

(4)

where 293.15 K is the room temperature, $L_{ti}$ is the tank insulation thickness.

This term is calculated individually for each layer based on specific tank insulation characteristics, including insulation thickness and exposure surface area. For intermediate layers, the area $A_{{wall}_i}$ represents the side area, while for the top and bottom layers, it incorporates both side and top/bottom base areas.

The term determining the relationship between water intake and heat balance is defined as:

$${\dot{Q}}_{hwc}=C_w\left(T_{w_{20}}-T_{cold}\right){\dot{m}}_{t,fraction},$$

(5)

where ${\dot{m}}_{t,fraction}$ is the fractional flow rate of the liquid.

The term ${\dot{Q}}_{hwc}$ uses the temperatures of the lower tank layer $T_{w_{20}}$ and incoming cold water $T_{cold}$ (15 °C) to evaluate the net impact of water intake on heat distribution, which predominantly affects the tank’s bottom two layers.

The liquid consumption rate for water intake is determined as follows:

$${\dot{m}}_{t,fraction}={\dot{m}}_t·fraction,$$

(6)

where ${\dot{m}}_t$ is the total consumption rate of hot water and the fraction is defined as:

$$fraction={\displaystyle \frac{\left(T_{use}-T_{cold}\right)}{\left(T_{w_0}-T_{cold}\right)}},$$

(7)

where $T_{use}$ is the required temperature (K) for using hot water (45 °C), $T_{w_0}$ is the upper layer of water from which the water is collected. In the current calculations, the fraction is considered constant and equal to 0.92.

Energy balance equations for Tank-2:

$$C_w{m}_{w_i}{\displaystyle \frac{d{T}_{w_i}}{dt}}={\dot{Q}}_{c,w_i}+{\dot{Q}}_{c,wf-w_j}-{\dot{Q}}_{{loss}_i}-{\dot{Q}}_{hwc},$$

(8)

where the main difference is the term ${\dot{Q}}_{c,{wf}_j-w_j}$ which represents the heat transfer rate between the heat exchanger and the liquid. Index $j$ denotes the layers of liquid where the heat exchanger is located. The heat exchanger, according to the tank configuration, is located from layer 9 to layer 18.

The heat transfer rate between the heat exchanger and the liquid in the tank is expressed as follows:

$${\dot{Q}}_{c,{wf}_j-w_j}=U_{wf-w_j}\left(T_{wf}-T_{w_j}\right)A_{{wall}_j},$$

(9)

where $U_{wf-w_j}$ is the overall heat transfer coefficient between the heat exchanger and the liquid in the tank, $T_{wf}$ is the temperature of the working fluid in the heat exchanger.

The overall heat transfer coefficient is calculated using the following equation:

$${\displaystyle \frac{1}{U_{wf-w_j}{A}_{{wall}_j}}}={\displaystyle \frac{1}{h_o{A}_o}}+R+{\displaystyle \frac{1}{h_i{A}_i}},$$

(10)

where $h_o$ is the heat transfer coefficient from the IHX outer surface, $h_i$ is the heat transfer coefficient from the IHX inner surface, $R$ is the resistance of the IHX pipe.

The heat transfer coefficient from the outer and inner surfaces are calculated as:

$$h_o={Nu}_n{\displaystyle \frac{k_w}{D_o}},$$

(11)

$$h_i={Nu}_{wf}{\displaystyle \frac{k_w}{D_i}},$$

(12)

The Nusselt number ${Nu}_n$ based on the pipe diameter is calculated using the correlation equation [22,23,24]:

$${Nu}_n=C_n·R{a}^n,$$

(13)

where coefficients $C_n$ and $n$ are constant values and represent a multiplier for natural convection, and an indicator of the degree of natural convection, respectively. $C_n$ and $n$ are taken to be $0.4$ and $0.25.$

The Rayleigh number is defined by the standard equation and is written as:

$$Ra={\displaystyle \frac{g{\beta }_w\left|T_{aver}-T_{w_i}\right|D_o^3}{{\nu }_w{\alpha }_w}},$$

(14)

where $D_o$ is the outer diameter of the IHX pipe, $T_{aver}$ is the average temperature between the liquid in the tank and the liquid in the IHX at a given calculation layer.

The Nusselt number from Equation (12) is calculated using two additional equations depending on whether the flow is laminar (Equation (15)) or turbulent (Equation (16)) [23,24] in the heat exchanger pipe. This equation and definition criterion were also taken from the TRNSYS 18 program [25]:

$${Nu}_{wf}={\left({\left({\displaystyle \frac{48}{11}}+{\displaystyle \frac{\raisebox{1ex}{$51$}\!\left/ \!\raisebox{-1ex}{$11$}\right.}{1+{\left(\frac{\raisebox{1ex}{$1342$}\!\left/ \!\raisebox{-1ex}{${Pr}_{wf}$}\right.}{H{E}^2}\right)}^2}}\right)}^3+1.816{\left({\displaystyle \frac{HE}{1+\frac{1.15}{{Pr}_{wf}}}}\right)}^{1.5}\right)}^{{\displaystyle \frac{1}{3}}},$$

(15)

$${Nu}_{wf}=0.023{Re}_{wf}^{0.85}{Pr}_{wf}^{0.4}{\left({\displaystyle \frac{D_i}{D_{coil}}}\right)}^{0.1},$$

(16)

In Equation (15), the coefficient $HE$ is determined by the following equation:

$$HE={\displaystyle \frac{{Re}_{wf}{\left(\frac{D_i}{D_{coil}}\right)}^{0.5}}{1+{\left(\frac{{Pitch}_{coil}}{\pi D_i}\right)}^2}},$$

(17)

where ${Pitch}_{coil}$ is the pitch of the heat exchanger coil.

The definition of the flow mode is calculated using the critical Reynolds number for round pipes:

$${Re}_{crit}=\mathrm{20,000}{\left({\displaystyle \frac{D_i}{D_{coil}}}\right)}^{0.32},$$

(18)

where $D_i$ is the inner diameter of the IHX tube and $D_{coil}$ is the diameter of the coiled tube heat exchanger. To check the mode, the current Reynolds number is calculated using the well-known equation:

$$Re={\displaystyle \frac{4\dot{m}}{\pi D_i{\mu }_{wf}}},$$

(19)

The Prandtl number is calculated using the equation:

$${Pr}_{wf}={\displaystyle \frac{C_{wf}{\mu }_{wf}}{k_{wf}}},$$

(20)

The following equation is used to calculate the charging efficiency:

$${\eta }_{ch}={\displaystyle \frac{\sum _t\dot{m}{C}_p\left(T_{in}-T_{out}\right)}{\sum _t\dot{m}{C}_p\left(T_{in}-T_{inital}\right)}},$$

(21)

where $T_{in}$ is the temperature of the water entering the tank or heat exchanger (45 °C), $T_{out}$ is the temperature of the last layer of water for both tanks, $T_{inital}$ is the initial temperature in the tank (20 °C).

For the discharge mode, the efficiency is calculated using the equation:

$${\eta }_{dis}={\displaystyle \frac{\sum _t\dot{m}{C}_p\left(T_{w_1}-T_{in}\right)}{\sum _t\dot{m}{C}_p\left(T_{inital}-T_{aver}\right)}},$$

(22)

where $T_{w_1}$ is the temperature of the upper layer of water that is taken, $T_{in}$ is the temperature of the cold water that is supplied for replacement (15 °C), $T_{aver}$ is the average temperature in the tank.

The Richardson number is calculated using the well-known equation:

$$Ri={\displaystyle \frac{g\beta \left(T_{w1}-T_{w20}\right)L_{wt}}{V^2}},$$

(23)

where $T_{w1}, T_{w20}$ are the temperatures of the upper and lower layers, respectively, $L_{wt}$ is the height of the tank, $V$ is the velocity of incoming hot water in Tank-1.

To estimate the error, the normalized root mean square error (RMSE) is used, which is calculated using the equation:

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}={\displaystyle \frac{1}{{T_w}^{max}-{T_w}^{min}}}\sqrt{{\displaystyle \frac{\sum _{j=1}^n\sum _{i=1}^t{\left(‖{T^{\times }}_{i,j}-T_{i,j}‖\right)}^2}{nt}} },$$

(24)

where ${T_w}^{max}, {{ T}_w}^{min}$ are the maximum and minimum temperatures in the tank for the entire calculation time, ${T^{\times }}_{i,j}$ is the measured temperature, $T_{i,j}$ is the numerically predicted temperature, $n$ is the number of layers, $t$ is the final calculation time. For the charging mode, the ${T_w}^{max}$ was set to 307.85 K, and the ${T_w}^{min}$ to 295.55 K. For the discharging mode, the corresponding values were 321.52 K and 295.65 K, respectively.

4. Calculation Method and Model Validation

To validate the numerical calculation algorithm, experimental data from Nash et al. [26] were used for comparison. The validation utilized a domestic hot WST with a capacity of 0.151 m³ (151 L), featuring an external height of 1.42 m, an internal height of 1.30 m, an inner diameter of 0.38 m, and an outer diameter of 0.51 m. An IHX with a diameter of 0.0254 m and a fluid capacity of 0.0075 m³ (7.5 L) is positioned at the bottom of the tank. Temperature measurements were obtained using thermocouples (TC) placed horizontally at the cold-water inlet and midway along the tank, and a custom TC wand at the hot water outlet. Experiments were conducted in a controlled environment at 23 °C. Data collected included fluid temperatures at the coil inlet and outlet, and volumetric flow rates for both the hot water discharge and the IHX water, with a sampling rate of approximately 0.15 s. Figure 2 presents a schematic representation of the test tank, as adapted from the work of Nash et al. [26].

The mathematical model outlined in Section 3, which is based on a system of ordinary differential equations (ODEs) with all associated coefficients, was numerically implemented in Python using the fourth-order Runge–Kutta method. The resulting simulations were compared with the experimental data reported by Nash et al. [26].

Figure 3 presents a comparison of tank temperatures between the experimental data and the numerical simulation for charging mode. The charging experiment was conducted over 120 min. During this phase, with a constant temperature of 50 °C supplied through the heat exchanger, temperatures were evaluated at two specific layers: the lowest and the highest. The figure demonstrates the proximity of these two layers to the experimental values. The normalized RMSE is 8.6% for the upper layer and 8.4% for the lower layer, with an overall normalized RMSE of 6.1%. These accuracy metrics are consistent with the performance of the model in the study by Nash et al. [26] for 20 layers in the tank. Thus, the results validate the reliability of the model during the charging mode.

Figure 4 compares temperatures obtained numerically with experimental data over 40 min during the discharging mode. In this mode, the heat exchanger was completely turned off, and only the water intake and addition of cold water at 24 °C were active. The experimental setup utilized seven sensors.

For comparison with the 20 layers used in the numerical model, data from layers 20, 18, 16, 13, 11, 6, and 3 were selected. The comparison, as shown in Figure 4, indicates good agreement between the experimental and numerical profiles. The overall normalized RMSE is 6.8%, which aligns with the expected accuracy for 20 layers, as discussed in Nash et al. [26]. The maximum deviation in the top two layers from the experimental results is attributed to diffusion error, arising from the approximation of the first-order time derivative. This supports the validity of the numerical model for the discharging mode.

5. Results and Discussions

This section presents a detailed numerical analysis of the charging and discharging behaviors of Tank-1 and Tank-2, examining how differences in heat exchanger configurations influence thermal stratification, charging efficiency, and discharging dynamics. During charging mode, variations in heating patterns are expected due to the unique heat exchanger arrangements. In the discharging mode, however, Tank-1 and Tank-2 will exhibit similar behavior, as both have identical cold-water feed and hot-water intake structures. For numerical modeling, the water in the tank is divided into 20 layers for both configurations. It is assumed that, at the start of the charging process, the entire tank’s water is at a temperature of 20 °C, while at the start of the discharging process, the water temperature is 45 °C.

5.1. Tank-1 During Charging Mode

Tank-1 is a system where hot water is supplied directly to heat the water inside the tank. The hot water is introduced at the top, while cold water is drawn from the bottom to promote efficient heating. The tank took 2.32 h, or approximately 2 h and 20 min, to reach full charge, with a measured charging efficiency of 94.57%. This configuration allows for efficient use of buoyancy-driven convection, as hot water naturally rises and cooler water is displaced downward. This process supports rapid thermal stratification, which is especially useful in applications requiring fast charging cycles, such as in solar-assisted water heating systems.

Figure 5 illustrates the temperature profile of the 20 layers over time. It is observed that the bottom layer is the last to experience a temperature increase due to the natural upward movement of water. This effect is attributed to buoyancy, as warmer liquid, being less dense, rises. Since hot water is introduced at the top of the tank, the temperature increases at the bottom layer more slowly. The analysis of Figure 5 suggests that the buoyant force has a greater impact than the flow rate of water supply and withdrawal. This is supported by the observation that the upper layers of the tank heat up significantly faster than the lower layers, despite the hot water being supplied and the water being withdrawn at equal rates.

This observation is supported by the Richardson number depicted in Figure 6, which illustrates its variation over time. During the initial 20 min, the Richardson number rises to its peak value due to the onset of the tank’s charging from its initial state. As the temperature increases, the buoyancy effect in the fluid also increases. The critical Richardson number of 400 indicates the ratio by which convective turbulence (buoyancy force) exceeds the mechanical effects (such as the speed of the hot water supply). As the temperature throughout the tank becomes more uniform, the Richardson number decreases because the reduction in temperature differences lowers the convective flows. This peak in the Richardson number demonstrates that the buoyancy forces initially dominate, establishing a strong thermal gradient that stabilizes over time. Understanding this dynamic can inform control strategies for systems with direct hot water input, optimizing the balance between flow rate and thermal stratification for different applications.

However, the presence of forced fluid movement allows for the observation of temperature distribution throughout the tank due to the artificial mixing of liquid layers. Figure 7 illustrates the variation in temperature distribution across layers over time. This figure demonstrates that the upper layers heat up significantly faster. Within 30 min, the hot water has risen to the top, resulting in the highest temperature in the upper layer. As a result of the forced circulation and intake, the temperature gradually equilibrates in the lower layers. This leads to the observable shift in the thermocline over time, as the lower layers adjust to the temperature of the upper layers. This dynamic is advantageous in cold climates, where maintaining a stable thermal gradient allows for consistent access to hot water in the upper layers, improving user comfort and energy efficiency.

5.2. Tank-2 During Charging Mode

The second type (Tank-2) is a tank utilizing an indirect charging method. Instead of directly supplying hot water to the tank, a heat exchanger is used. This heat exchanger is a coil located within the tank at a height ranging from 0.15 m to 0.6 m. The charging process lasted 3.78 h, or approximately 3 h and 47 min, with a charging efficiency of 85.87%.

Figure 8 illustrates the temperature evolution of the layers over time. The lower layer heats up much more slowly, which is attributed to the heat exchanger not covering all the lower layers of the tank. This results in the formation of a pronounced thermocline, which remains stable throughout the charging process. The thermocline is a consequence of buoyancy within the fluid, as there is no forced mixing due to the absence of direct hot water flow. Therefore, the Richardson number is not applicable in this context, as the fluid movement within the tank is driven solely by natural convection (buoyancy forces).

Figure 8 also shows that the highest temperatures are observed in the layers where the heat exchanger is present, followed by the layers located above the heat exchanger. The lowest temperatures are found in the bottom layer, which is not influenced by the heat exchanger. This occurs because the layers in contact with the heat exchanger heat up more rapidly, and the heat gradually rises to the upper layers due to buoyancy forces. The temperature in the layers where the heat exchanger is located remains relatively higher during charging, as these layers continue to receive heat directly from the exchanger.

Figure 9 shows the temperature variation along the height of the tank. Higher temperatures are observed in the regions where the heat exchanger is located. Although the rapid heating in these regions might seem non-physical, the final calculations of this model align with the experimental results reported by Nash et al. [26]. The distribution shows that the upper layers are slowly heated by convective flows originating from the layers in contact with the heat exchanger, leading to a gradual temperature equalization across the tank.

The thermocline is clearly established just above the heat exchanger, where the high temperatures create a distinct separation between the upper and lower layers. Despite this, convection within the layers results in gradual, chaotic mixing, causing the lower layers to heat up slowly over time. This gradual heating of the lowest layer is both physically accurate and experimentally validated.

5.3. Tank-1 and Tank-2 During Discharging Mode

The discharging process for both tanks is identical. Hot water is drawn from the top while cold water is introduced at the bottom at a rate of 330 kg/h. The discharge lasted 0.55 h, or approximately 33 min, with a discharging efficiency of 13.63%.

Figure 10 illustrates the temperature variation across the layers over time during discharge. As cold water is introduced at the bottom, the lower layer experiences a rapid decrease in temperature. Over time, the upper layers also cool down. This figure shows the thermocline gradually rising, now driven by the influx of cold water rather than hot. As the volume of hot water diminishes, the cold temperature progressively spreads upward, cooling the tank from the bottom to the top.

Figure 11, which depicts the temperature variation by height, more accurately illustrates the shift of the thermocline. The distribution shown in Figure 11 aligns well with the experimental data. Additionally, a similarity can be observed between the temperature distribution during charging in Figure 7 for Tank-1 and during discharging in Figure 11. In Tank-1, charging is achieved through the direct addition of hot water, which is analogous to the introduction of cold water during the discharging process. Consequently, the dynamics of temperature change by height are inversely similar in these two cases, with charging resembling the reverse process of discharging, where hot water is drawn from the tank.

Overall, these findings highlight the distinct advantages of each configuration based on specific use cases. Tank-1’s rapid charging and responsive discharging make it ideal for high-demand applications such as residential hot water systems in colder climates. Meanwhile, Tank-2’s stable thermal profile supports applications that benefit from prolonged heat retention, such as in storage systems used to supplement heating during evening hours. Future studies could further investigate how variable inlet temperatures and flow rates affect these stratification patterns, providing insights into optimizing TES systems for a wide range of climatic and operational conditions.

6. Conclusions

This study conducted a comparative numerical analysis of two thermal energy storage (TES) tank configurations: one with an external heat exchanger (Tank-1) and the other with an internal heat exchanger (Tank-2). The results demonstrate that Tank-1 achieves a higher charging efficiency of 94.57%, reaching full charge in approximately 2 h and 20 min, whereas Tank-2 requires about 3 h and 47 min with a charging efficiency of 85.87%. Both tanks exhibit similar performance during the discharging process, achieving a discharging efficiency of 13.63% over 33 min.

The analysis reveals that the external heat exchanger configuration in Tank-1 promotes effective thermal stratification and achieves a more uniform temperature distribution, influenced by buoyancy effects, as shown by Richardson number analysis. This design is particularly suitable for applications requiring rapid heating and minimal heat loss, making it beneficial in colder climates or systems with intermittent demand. Conversely, the internal heat exchanger in Tank-2 supports a stable but less efficient stratification, which is advantageous for applications needing a gradual and sustained release of heat, such as prolonged heating in moderate climates.

These findings highlight the potential for Tank-1’s external heat exchanger design in environments where fast and efficient thermal energy storage is essential, while Tank-2 could be more beneficial in contexts where steady heat retention is prioritized. The adaptive-grid model and the detailed energy balance equations implemented in this study have shown reliable predictive accuracy, with an RMSE of 6.1% during charging and 6.8% during discharging, validating the model against experimental data.

Future work could focus on integrating these tank configurations with specific types of heating systems, such as heat pumps, solar collectors, or latent heat thermal energy storage systems using phase change materials (PCMs). Integrating the studied tanks with renewable energy technologies and energy efficiency solutions could help smooth the intermittent nature of these systems and optimize supply-demand curves. Additionally, analyzing the impact of variable flow rates, inlet temperatures, and different heat transfer fluids could further enhance understanding of TES tank performance under diverse environmental and operational conditions. Future research could also explore lifecycle and economic assessments to evaluate the cost-benefit ratios of each configuration for specific applications, ensuring that TES systems are both efficient and economically viable for widespread deployment.