# A hybrid model for packed bed thermal energy storage system

Shri Balaji Padmanabhan, Mohamed Tahar Mabrouk, Bruno Lacarrière

► To cite this version:

Shri Balaji Padmanabhan, Mohamed Tahar Mabrouk, Bruno Lacarrière. A hybrid model for packed bed thermal energy storage system. Journal of Energy Storage, 2024, 98, pp.113068. 10.1016/j.est.2024.113068. hal-04746406

HAL Id: hal-04746406

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

Submitted on 21 Oct 2024

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
<figure><img src="image_2.png" /></figure>

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

# A hybrid model for packed bed thermal energy storage system

Shri Balaji Padmanabhan *, Mohamed Tahar Mabrouk *, Bruno Lacarrière

IMT Atlantique, Department of Energy Systems and Environment, GEPEA UMR CNRS, 6144, F-44307, Nantes, France

## ARTICLE INFO

Dataset link: https://gitlab.imt-atlantique.fr/hybrid_modeling/heat-storage

**Keywords:**

- Thermal Energy Storage
- Hybrid modeling
- Numerical model
- Machine learning
- Deep Neural Networks

## ABSTRACT

Physical models for Packed Bed Thermal Energy Storage (PBTES) system plays a crucial role in predicting its dynamic behavior and long-term performance under various operational conditions. Nevertheless, the prevailing physical models struggles to balance between computational efficiency and accuracy. This paper introduces a novel hybrid model for PBTES that combines the strengths of numerical modeling and machine learning, aiming to address this challenge. The proposed model includes a low-precision yet faster-to-solve linearized version of the state-of-the-art two-phase model, integrated with a machine learning module. In the hybrid model, the linear two phase model is solved numerically using finite volume method on a coarse mesh. Then, the ML module takes this low precision and linear solution, and maps it into a high precision non-linear solution corresponding to fine mesh. The proposed hybrid model is highly robust and generalized, and has been found to be effective in handling cases even beyond the training range of the machine learning module. During the validation process, the hybrid two-phase model delivered the non-linear solutions nearly 350 times faster than the traditional non-linear two-phase model. It demonstrated good accuracy, achieving an overall error percentage just 0.16 and a R² score of 0.99. The findings highlight the potential of this hybrid modeling approach for fast and accurate simulations of PBTES, suggesting its applicability in complex applications such as model control optimization and long period simulations.

## 1. Introduction

Thermal Energy Storage (TES) is pivotal in the energy sector. It contributes significantly towards improved energy efficiency and reliability of energy systems and reduces their environmental impacts. Specifically, TES enables increased integration of renewable energy sources by mitigating their intermittency. Also, it serves to reduce the effects of fluctuations in the demand on the production, and functions as a thermal buffer. It operates by storing thermal energy during periods of low demand or excess generation in the form of sensible heat or latent heat. The stored energy can then be released as needed, which can help to balance the demand and supply of thermal energy.

In concentrated solar power plants, TES can be used to store solar energy surplus effectively and ensures continuous power generation by using the stored energy during periods of limited solar availability and also helps in the development of solar thermal systems [1,2]. In adiabatic compressed air energy storage, integration of TES increases the storage efficiency and ensures a stable and on-demand thermal output [3]. In the industrial context, TES improves energy efficiency by capturing, storing, and recovering the waste heat produced during the production process [4,5]. In district heating networks, TES can store excess heat during periods of low demand and helps in load balancing

by releasing it during peak demands. Also, it helps in better integration of renewable energy sources into the district heating system.

One of the most cost-effective variants of TES is Packed Bed Thermal Energy Storage (PBTES) [6]. It consists of a container filled with a solid filler material, through which a heat transfer fluid is circulated. It can store energy in the form of sensible or latent heat. Sensible heat storage is known for its superior safety and cost-effectiveness and, it is mostly preferred for both low and high-temperature heat storage applications. Conversely, latent heat storage, which involves phase change materials, typically more expensive and is more suitable for applications necessitating higher storage densities and low to medium temperature levels [7]. PBTES systems are noted for their versatility and finds a wide range of applications across the energy sector.

The dynamic behavior and performance of PBTES are governed by various complex heat transfer interactions. Therefore, numerical models are employed to facilitate the modeling, control, and design optimization of PBTES by predicting its dynamic behavior. In the literature, several numerical models have been developed to model the dynamic behavior of PBTES. There are three major numerical models for PBTES, which includes the single phase model [8], the two phase model [9], the three phase model [10]. The two phase model

<footnote>* Corresponding author.
E-mail addresses: shri-balaji.padmanabhan@imt-atlantique.fr (S.B. Padmanabhan), mohamed-tahar.mabrouk@imt-atlantique.fr (M.T. Mabrouk).</footnote>
Nomenclature

A Specific surface area (m²)

Cₚ Specific heat capacity (JK⁻¹kg⁻¹)

D Diameter (m)

k Thermal conductivity (Wm⁻¹K⁻¹)

T Temperature (K)

U Fluid velocity (ms⁻¹)

t Time coordinate (s)

z Spatial coordinate (m)

L Tank length (m)

M Mass (kg)

V Volume (m³)

ṁ Mass flow rate (kgs⁻¹)

Pe Peclet number (-)

Sᵢ Time scaling factor (s)

Bi Biot number (-)

Re Reynolds number (-)

h Heat transfer coefficient (Wm⁻²K⁻¹)

hᵥ Volumetric heat transfer coefficient (Wm⁻³K⁻¹)

CV Control Volume

H Volumetric heat transfer coefficient between tank and ambient (Wm⁻³K⁻¹)

Greek letters

ε Void fraction (-)

μ Dynamic viscosity (kgm⁻¹s⁻¹)

ρ Density (kgm⁻³)

θ Dimensionless temperature (-)

ξ Dimensionless spatial coordinate (-)

τ Dimensionless time (-)

Subscripts

0 Initial

int Internal

f Fluid

s Solid

sf Fluid and solid shared property

in Inlet

env Ambient

max Maximum

is the most popular option and is widely adopted in previous studies on numerical simulation of PBTES [7,9,11,12]. The two phase model generally provides accurate approximation of the thermal performance due to its consideration of temperature gradients and heat transfer dynamics between the thermal fluid and the solid filler material. But the down side is that it is computationally expensive and hence, it is time-consuming to solve this model [13]. Consequently, they are often ill-suited for optimal model control, long time period simulations and other real-time applications. To overcome the computational burden of two phase model, a perturbation model [14] has been developed. There have also been efforts to develop analytical [11,15] and semi-analytical [16] solutions for the perturbation model. While the perturbation model and its analytical forms offers a simplified and computationally efficient alternative to the two phase model, it has limitations in handling variable inputs and non-linear properties due to its underlying assumptions and approximations leading to a less accurate model.

As Artificial Intelligence techniques continue to develop rapidly, the prospect of using Machine Learning (ML) in the modeling of PBTES is becoming more interesting and promising. ML models are computationally efficient and provide faster predictions compared to traditional numerical models. They are robust and accurately predicts key performances by effectively managing non-linearity and complex interactions in PBTES [17]. For the purpose of optimization, data driven ML models for thermal storage have been developed with neural networks [18–20], ensemble learning models [21] and neural networks with Group Method of Data Handling (GMDH) [22], which provides faster prediction of thermal performances and facilitates faster optimization process. Apart from these studies, ML models for PBTES are not extensively studied in the literature. Even though, ML models provide faster predictions with good accuracy, they are heavily reliant on the quantity and quality of the training data and prone to overfitting. Also they have weak generalization which may lead to poor performance outside their training range, resulting in a reduction of confidence in their predictions [17]. Furthermore, the major disadvantage of these models is that they often do not incorporate underlying physical principles, making their predictions purely data-driven and result in patterns that adhere to the training data but defy the known laws of physics outside that range.

One feasible approach is the implementation of hybrid models which combines the benefits of physics-based and machine learning

models. These models constitute a blend of numerical and machine learning methodologies and they comply with physical laws while simultaneously exhibit high accuracy and speed. Hybrid modeling has already made a significant impact across a wide range of scientific disciplines such as tomography image processing [23], environmental science [24], water systems modeling [25], energy processes [26], robotics [27], chemical reactions modeling [28], and groundwater modeling [29]. This class of methods harness the strengths of both numerical models and machine learning models, and they have improved model transparency, explainability, and analytic capabilities at reduced computational cost.

This research primarily aims to improve the state-of-the-art two-phase model by making it faster and easier to solve, while maintaining accuracy through a hybrid modeling approach. In this paper, a hybrid Two Phase model for sensible heat version of PBTES is proposed. The hybrid version of the two phase model presented here combines a linearized version of two phase model solved using a coarse mesh with a machine learning module to generate high precision non-linear solutions corresponding to fine mesh. The usage of coarse mesh does not compromise the accuracy of the hybrid model, and during validation, the hybrid model demonstrated excellent computational performance with high accuracy.

This manuscript is organized as following: Section 2 describes the state of art mathematical model of the PBTES, details the numerical method used to solve it and highlights the difference between the “linear” and “non-linear” variants of the model in terms of computational time and accuracy. Section 3 describes the novel hybrid model of the PBTES. Section 4 presents the model’s performance on various validation sequences, including those covering operational conditions beyond the training set, and discusses the results. Section 5 concludes the paper.

## 2. Mathematical model and numerical methodology

The storage tank, as illustrated in Fig. 1, consists of a vertical cylinder filled with solid filler material and a heat transfer fluid which flows through the filler material. The top and bottom of the tank have distributors with ports, which ensures that the incoming flow is uniformly distributed over the cross-sectional area of the tank. There are several options for the solid filler material such as rocks, pebbles, metals and ceramics. Mostly, they are recycled materials from waste.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Schematic diagram of packed bed thermal energy storage.</figcaption></figure>

For heat transfer fluid, the most popular options are thermal oils, steam and water. The operation of the PBTES comprises two processes namely charging and discharging. During the charging process, hot fluid enters the tank through the top port and diffuses through the filler material, transferring heat from the fluid to the filler material, and exits from the bottom port. Conversely, during the discharging process, the flow direction is reversed, the cold fluid enters from the bottom port, and as it passes through the hot filler material, heat is transferred from the filler material to the fluid and hot fluid exits from the top.

## 2.1. Numerical model

In this study, state of the art Two-Phase (TP) model for PBTES is utilized. This model has been validated against experimental results in previous studies [9,30], demonstrating its robustness and reliability. The following assumptions were made in the TP model,

(i) The inlet fluid velocity is assumed to be uniform along the radial axis, so the distributors are not modeled.

(ii) The fluid is incompressible, and the momentum equations are not modeled.

(iii) The solid is represented as a continuous, homogeneous, and isotropic granular medium.

(iv) Thermal diffusion is considered only along the tank's axial direction, resulting in one-dimensional energy equations.

With the above assumptions, the one-dimensional energy equations for fluid and solid as provided in [9], are given by Eqs. (1) and (2) respectively.

$$ \epsilon (\rho C_p)_f \frac{\partial T_f}{\partial t} + \epsilon (\rho C_p)_f U \frac{\partial T_f}{\partial z} = \frac{\partial}{\partial z} \left( k_f \frac{\partial T_f}{\partial z} \right) + h_v (T_s - T_f) - H (T_f - T_{env}) \quad (1) $$

$$ (1 - \epsilon) (\rho C_p)_s \frac{\partial T_s}{\partial t} = k_s \frac{\partial^2 T_s}{\partial z^2} - h_v (T_s - T_f) \quad (2) $$

Boundary and initial conditions of Eqs. (1) and (2) are given as following,

$$ \begin{aligned} T_f(t, z=0) &= T_{f,\text{in}}; \quad \left(\frac{\partial T_f}{\partial z}\right)_{z=L} = 0 \\ \left(\frac{\partial T_s}{\partial z}\right)_{z=0} &= 0; \quad \left(\frac{\partial T_s}{\partial z}\right)_{z=L} = 0 \\ T_f(t=0, z) &= T_{f,0}(z); \quad T_s(t=0, z) = T_{s,0}(z) \end{aligned} $$

The above boundary conditions corresponds to the charging process and for the discharging process the fluid boundary conditions are reversed. The velocity of the fluid U is given by Eq. (3). The volumetric heat transfer coefficient $h_v$ between fluid and solid is determined using the correlation provided by [31], which is given by Eq. (4). Furthermore, to account for the heat capacity contribution of the storage wall in the model, solid density is modified by integrating the wall's heat capacity into the calculation. The modified solid density $\hat{\rho}_s$ is given by Eq. (5). The addition of wall contribution leads to a slower numerical evolution due to add-on thermal inertia, but makes TP model's solution align better with experimental results [9].

$$ U = \frac{\dot{m}}{\rho_f \pi \frac{D_v^2}{4} \epsilon} \quad (3) $$

$$ h_v = \frac{6(1 - \epsilon)k_f(2 + 1.1 R e_f^{0.6} P r_f^{0.33})}{D_{\text{rock}}^2} \quad (4) $$

$$ \hat{\rho}_s = \rho_s + \frac{M_{\text{wall}} C_p \omega_{\text{wall}}}{(1 - \epsilon) C_p s V_{\text{int}}} \quad (5) $$

To enhance the generality and applicability of the study, the TP energy equations are converted into a dimensionless form and then analyzed. This transformation allows the results to be scaled for different applications, ensuring the model's scalability for both small and large-scale systems. The TP energy equations are transformed into a non-dimensional form using following transformations,

$$ \theta = \frac{T - T_{\text{min}}}{T_{\text{max}} - T_{\text{min}}}; \quad \xi = \frac{z}{L}; \quad \tau = \frac{t}{S_t}; \quad S_t = \frac{L}{\epsilon U_{\text{max}}} $$

The time t and spatial coordinate z are scaled with a time scaling factor $S_t$ and the length of the storage tank L respectively. Furthermore, the temperatures of the fluid and solid are normalized within the cold ($T_{\text{min}}$) and hot ($T_{\text{max}}$) operating temperatures of the thermal storage. In this study, $T_{\text{min}} = 100$ °C and $T_{\text{max}} = 250$ °C are used as operating temperature range of PBTES. The non-dimensional form of the TP energy equations are given by Eqs. (6) and (7) respectively.

$$ \frac{\partial \theta_f}{\partial \tau} + \text{Pe} \frac{\partial \theta_f}{\partial \xi} = \alpha_f \frac{\partial}{\partial \xi} \left( k_f \frac{\partial \theta_f}{\partial \xi} \right) + \beta_f (\theta_s - \theta_f) - \gamma_f (\theta_f - \theta_{\text{env}}) \quad (6) $$

$$ \frac{\partial \theta_s}{\partial \tau} = \alpha_s \frac{\partial^2 \theta_s}{\partial \xi^2} - \beta_s (\theta_s - \theta_f) \quad (7) $$

The non dimensionless parameters used in the transformation are given as following,

$$ \text{Pe} = \frac{U S_t}{L}; \quad \alpha_f = \frac{S_t}{\epsilon(\rho C_p)_f (L)^2}; \quad \beta_f = \frac{h_v S_t}{\epsilon(\rho C_p)_f}; \quad \gamma_f = \frac{H S_t}{\epsilon(\rho C_p)_f} $$

$$ \alpha_s = \frac{k_s S_t}{(1 - \epsilon)(\rho C_p)_s (L)^2}; \quad \beta_s = \frac{h_v S_t}{(1 - \epsilon)(\rho C_p)_s} $$

The dimensionless parameters $Pe$ and $\beta$ represents the flow characteristics, and $\alpha$ represents the material properties of the system. The relative importance of heat transfer between the fluid and the ambient environment compared to the storage capacity of the fluid is represented by $\gamma$.

## 2.2. Numerical method

The TP model can be solved numerically using the finite volume method [30,32–34], the finite difference method [9], or method of characteristics [12]. Among these methods, Finite Volume (FV) method is a powerful numerical technique. Its primary strength lies in the conservation of fluxes across the control volume (CV)'s boundaries, ensuring energy conservation. It can handle discontinuities and thermal gradients with ease making it robust and accurate. This explains its popularity in the literature. Therefore, in this study, the transient TP model is solved using the finite volume method on a uniform grid. In FV method, the storage domain is divided into n number of CVs and the dimensionless equations (6) and (7) are integrated over each CV, resulting in 2 * n number of Ordinary Differential Equations (ODEs). Then, ODEs are solved using ODE solvers for time integration.
<table><caption>**Table 1**
Thermo-physical properties of fluid and solid.</caption><thead><tr><th>Properties</th><th>Linear model</th><th>Non-linear model</th></tr></thead><tbody><tr><td>ρ<sub>f</sub></td><td>845</td><td>-0.000321(T<sub>f</sub><sup>2</sup>) - 0.614254(T<sub>f</sub>) + 1020.62</td></tr><tr><td>μ<sub>f</sub></td><td>0.00057</td><td>exp(-2.048 log(T<sub>f</sub>) + 10.773) × 1.0e<sup>-03</sup></td></tr><tr><td>k<sub>f</sub></td><td>0.1</td><td>-1.5e<sup>-07</sup>(T<sub>f</sub><sup>2</sup>) - 3.3e<sup>-05</sup>(T<sub>f</sub>) + 0.118294</td></tr><tr><td>Cp<sub>f</sub></td><td>2380</td><td>8.970785e<sup>-4</sup>(T<sub>f</sub><sup>2</sup>) + 3.313(T<sub>f</sub>) + 1496.005</td></tr><tr><td>ρ<sub>s</sub></td><td>2595</td><td>2595</td></tr><tr><td>k<sub>s</sub></td><td>5.5</td><td>5.5</td></tr><tr><td>Cp<sub>s</sub></td><td>900</td><td>0.8841(T<sub>s</sub>) + 795.9</td></tr></tbody></table>

## 2.3. Thermo-physical properties of fluid and solid

In the TP model, the thermo-physical properties of the fluid and the solid are generally dependent on their respective temperatures, making the TP model non-linear. The variations of thermo-physical properties with respect to temperature are generally expressed as polynomial equations. These temperature-dependent properties provides more realistic and accurate approximation of the behavior of the model by capturing the dynamic changes more effectively and making the results align better with experimental results [9,30]. Conversely, the thermo-physical properties of the model can be considered as temperature independent that is constants, leading to a linear TP model. This can be done by finding the average values of the thermo-physical properties in the range of the working temperature. Furthermore, in the linear model, since the thermo-physical properties of the fluid are not temperature dependent, it leads to a linear volumetric heat transfer coefficient. The linearized model is computationally efficient and easier to solve. However, this simplification overlooks critical variations in thermo-physical properties due to changes in temperature, resulting in discrepancies in predicting temperature profiles within the storage. This less computationally intensive linear TP model can be useful for preliminary studies, intense optimization, and simulation scenarios where high precision is not critical.

In this study a storage tank of length $L = 3$ m and diameter $D_t = 1$ m is considered. A mixture of silica gravel with diameter $D_{rock} = 0.03$ m and silica sand with diameter $D_{sand} = 0.003$ m in a mass proportion of 80% and 20% respectively is chosen as the solid filler material with a void fraction of 0.27. Therminol 88 that supports a maximum temperature of 350 °C is selected as working fluid. The thermo-physical properties of the solid filler material (gravel+sand) and the fluid (Therminol 88) used in both the linear and non-linear models in this study are listed in Table 1. The tank dimensions, filler material characteristics, the choice of the working fluid, and thermo-physical properties are in accordance with [9,30].

## 2.4. Grid convergence study

To determine the appropriate number of CVs needed for solving the non-linear TP model using the finite volume (FV) method, a grid convergence study is conducted. This study ensures the accuracy and stability of the numerical solution, verifying that the discretized model adequately captures the complex behavior of the TP model. The non-linear TP model is solved using FV method with different number of CVs. The fluid and solid solutions obtained for different number of CVs are compared against the respective reference solutions. The reference solution is obtained by solving the model using 2000 CVs, and it is considered grid-independent and highly accurate. For the convergence study, a simulation which comprise 4 consecutive charging and discharging process is performed. Mean Squared Error (MSE) is used as evaluating metrics. The variation of MSE and computation time for different number of CVs is given in Fig. 2. The grid convergence plot shows that beyond 1000 CVs, further refinement yields only minor improvements in accuracy, while significantly increasing computational time. Therefore, 1000 CVs have been chosen as the optimal number.

<table><caption>**Table 2**
Comparison of computation time - linear vs non-linear TP model.</caption><thead><tr><th>Model</th><th>Computation time [s]</th></tr></thead><tbody><tr><td>Linear</td><td>0.46</td></tr><tr><td>Non-linear</td><td>17.38</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Grid convergence: Variation of MSE and computation time over different number of CVs.</figcaption></figure>

## 2.5. Comparative study between linear and non-linear TP model

To showcase the effect of linearization of the non-linear TP model, a comparative analysis is conducted between the numerical solutions of the linear and non-linear models, providing a comprehensive understanding of the dynamics of their respective solutions under these defined conditions. In this study, the discretized FV ODEs for linear and non-linear TP model are solved in Julia programming language using the DifferentialEquations.jl [35] library's ODE solver.

Both linear and non-linear TP models are solved for the time varying inlet fluid temperature $\theta_{f,in}$ and inlet mass flow rate $\dot{m}_{in}$. The variation of $\theta_{f,in}$ is gradual and smooth over the time, where as $\dot{m}_{in}$ is varied in a step like manner over the time. The comparative study involves solving the models for two consecutive pairs of charging and discharging processes, and each processes is simulated for a dimensionless time $\tau = 15$. The time varying inputs for which the simulation performed is as depicted in Fig. 3. The solid temperature distributions solution along the storage length at different dimensionless times $\tau = 4, 8, 12$ are shown in Figs. 4(a) and 4(b) respectively. The fluid outlet temperature solution during each discharging process is shown in Figs. 5(a) and 5(b). Additionally, Fig. 6 illustrates the error between spatial solutions at each time step between non-linear and linear solutions for the entire time period of the simulation. This is done by calculating MSE between the corresponding non-linear and linear temperature distribution solutions in °C (full dimension) along the length of storage for each time step.

From the error distribution (Fig. 6), it is evident that a significant difference exists between the temperature solutions obtained by solving the linear and non-linear TP models. This discrepancy may keep-on continuing as the charging and discharging cycles progress. For the Linear TP model, the comparative study reveals that linearization does not adequately capture the proper dynamics of both the fluid and solid. From Tables 2 and 3, though the linear TP model is computationally efficient and capable of producing solutions quickly, it tends to have poor accuracy. In contrast, the non-linear TP model is almost 38 times slower than linear TP model but provides more precise approximation.

## 3. Hybrid two phase model

To bridge the accuracy gap between the computationally efficient linear TP model and more accurate but expensive non-linear model, a
<figure><img src="image_6.png" /><figcaption>Fig. 3. Input sequences for the comparative study: (a) Inlet fluid temperature, (b) Inlet fluid mass flow rate.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. Comparison of solid temperature distribution at three different instance during the charging processes ($\tau = 4, 8, 12$): (a) 1st Charging Process, (b) 2nd Charging Process.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Fluid outlet temperature during the discharging processes over the respective time periods: (a) 1st discharging Process and (b) 2nd discharging Process.</figcaption></figure>

hybrid TP model is proposed. It combines the linear TP model with a Machine Learning (ML) module to output high precision non-linear solution. Furthermore, in the hybrid model, the linear TP model is solved

using a smaller number of CVs ($n_{hyb}$) than the 1000 CVs established in the FV model convergence study. This initial step produces a coarse linear solution. The ML module takes this coarse and low precision
<table><caption>**Table 3**
Overall MSE between linear and non-linear TP model's solution.</caption><thead><tr><td>Phase</td><td>Dimension less [–]</td><td>Full dimension [°C]</td></tr></thead><tbody><tr><td>Fluid</td><td>2.00e-03</td><td>45.03</td></tr><tr><td>Solid</td><td>2.41e-03</td><td>54.22</td></tr></tbody></table>

linear solution, and maps it into a high precision non-linear solution that matches the numerical resolution of a solution corresponding to very fine mesh (1000 CVs). As this hybrid modeling approach primarily focus on learning the mapping from linearity to non-linearity, this hybrid modeling approach is easily extendable to unstructured and higher-dimensional meshes, as well as to other numerical methods.

### 3.1. Machine learning module

As previously mentioned, the ML module principally learns the mapping from coarse linear solution to fine non-linear solution. The use of neural networks to approximate fine grid solutions from coarse grid solutions has been explored in the literature. Neural networks have been employed to generate high-resolution mesh solutions from coarse or interpolated coarse solutions, demonstrating their potential for enhancing resolution in computational models [36,37]. In the context of mapping non-linear solution using neural networks, some studies in the literature have explored this by using the neural network to learn an implicit transformation that makes the non-linear system behave linearly in a new space [38], and transformation of large linear systems into small non-linear system using neural networks [39]. However, these approaches majorly differs from the methodology employed in this hybrid TP model, where non-linear model is explicitly linearized and ML module effectively performs mapping from linear solution to non-linear solution while also resolving the coarse to fine mesh solution at same time.

The ML module consists of a feature extractor, a Deep Neural Network (DNN) and a feature reconstructor. The feature extractor and reconstructor, plays a huge role by aiding the mapping process done by DNN. The feature extractor processes the coarse linear solution for each time step and extracts the dominant linear features that characterize the coarse solution. Since the solution is obtained by solving the linear model for fewer number of CVs, it may exhibits numerical oscillations or instabilities. This phenomenon occurs due to insufficient spatial resolution, which may fail to accurately capture the gradients and dynamics of the underlying physical problem. If this solution is directly given to DNN, the presence of noise and oscillations can make the training process unstable, causing convergence issues. Also, it prevents the neural network from learning the essential characteristics, leading to poor generalized and over-fitted model. These makes the feature extraction process significant. Then, DNN takes those extracted linear features (solid and fluid) along with primary inputs of linear model and maps into corresponding non-linear features (solid and fluid). Finally, the feature reconstructor, which reconstructs high resolution non-linear spatial solutions. Here, the primary computational cost lies mainly in solving the linear TP model. This hybrid model effectively streamlines the process, balancing the computational efficiency of linear TP model with the precision of non-linear TP model. Fig. 7 provides a visual representation of the hybrid TP model with detailed workflow.

### 3.2. Comprehensive design specifications of the ML module

In the ML module, Principal Components Analysis (PCA) [40,41] is utilized behind the feature extractor and reconstructor. PCA transforms the original values into a new set of uncorrelated values called principal components. Each principal component is a linear combination of the original values and captures as much variance as possible, with the first component capturing the most variance. The feature extractor takes in the linear temperature distribution solution corresponding to

<table><caption>**Table 4**
Architecture of the DNN.</caption><tr><td>No. of layers</td><td>6</td></tr><tr><td>No. of neurons per layer</td><td>100</td></tr><tr><td>Activation function</td><td>Leaky Rectified Linear Unit</td></tr></table>

$n_{hyb}$ CVs for each time step. Then it applies PCA transformation to this linear solution and extracts the features which is of reduced size $N_{pca}$. The choice of $N_{pca}$ is flexible and depends on $n_{hyb}$. It should be of an appropriate size to avoid noise and prevent carrying over oscillations from the coarse solution. Choosing a very small $N_{pca}$ can make the training process faster but it can lead to poor performance of ML module and vice-versa. This feature extraction process helps the DNN in effectively learning the map. The second component is DNN. For this context of mapping from coarse linear to refined non-linear solution, a simple DNN is found to be sufficient and performs well. The specifications of the DNN are very basic without any dropout and special regularization, and the details about architecture used in this study are outlined in Table 4. The last component in the ML module is the feature reconstructor. It takes the mapped non-linear features of size $N_{pca}$ and applies the inverse transformation of the PCA, and reconstructs the refined non-linear solutions corresponding to 1000 CVs.

### 3.3. Data generation and training of the ML module

The key input variables for the TP model are the working fluid's inlet mass flow rate and inlet temperature. The data generation process for this study involves a comprehensive simulations of both charging and discharging processes in linear and non-linear TP models. The linear TP model is solved with $n_{hyb}$ number of CVs (coarse mesh) and non-linear TP model is solved with 1000 number of CVs (fine mesh). The choice of $n_{hyb}$ will be discussed later. For the charging process, the inlet fluid temperature is varied between 0.6 to 1.0 in dimensionless value, which translates into 190 °C and 250 °C and the inlet mass flow rate is varied between 0.45 to 0.65 kg/s. Similarly, in the discharging process, the inlet fluid temperature is varied between 0.0 to 0.2 in dimensionless value, which translates into 100 °C and 130 °C and the inlet mass flow rate range is kept same as charging. For both processes, the inputs are gradually and smoothly varied in a completely random way over the time. Each charging and discharging process made to last for a dimensionless time duration of $\tau = 10$ (around 4 h in actual time). It is crucial to maintain the same time step value $\Delta t$ while solving the numerical models to generate training data, as this helps the DNN capture patterns between different input samples. However, during the testing and deployment phases, it is not necessary to maintain the same $\Delta t$, and adaptive time-stepping schemes can be used.

To generate data for training of ML module, 10 sequences of simulations are performed. In each sequence, 4 consecutive pairs of charging and discharging processes are continuously simulated. For every sequence, the initial charging process commenced with zero initial condition for both the fluid and the solid phase. Then, after each charging or discharging process, the final solution from the non-linear model of the preceding process was applied as the initial condition to both the linear and non-linear models for the current process. This is mainly done to mimic the workflow of the hybrid TP model in the prediction and deployment stage, where non linear mapping is performed after each process, and the final state of the mapped non-linear solution of the preceding process was applied as the initial condition for the current process. If the training data is not generated by following this technique, the difference between the linear and non-linear model will keep on progressing in a random manner without any pattern and makes the learning of mapping very difficult. Overall, this training data generation technique is crucial in minimizing error propagation in the linear model after each process and improving the machine learning
<figure><img src="image_9.png" /><figcaption>Fig. 6. Time step wise MSE between linear and non-linear solutions in full dimension.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. Schematic workflow of the hybrid TP model.</figcaption></figure>

model's ability to effectively learn the non-linearity during training phase.

The training of ML module is as follows: Separate PCA models are built for fluid and solid linear solutions, as well as for fluid and solid non-linear solutions. Thus, in total four PCA models are trained. The feature extractor uses the transformation part of the PCA models trained on the fluid and solid linear solutions to reduce the dimensions from $n_{hyb}$ to $N_{pca}$. Conversely, the reconstructor uses the inverse transformation component of the PCA models for fluid and solid non-linear solutions, transforming the data back from $N_{pca}$ to the original 1000 CVs. This approach ensures that the feature extraction and reconstruction processes are specifically tailored to the characteristics of respective solutions. The inputs for training the DNN are prepared by combining the extracted fluid and solid features with other inputs, and stacking them for each time step. The DNN is then trained to generate the corresponding fluid and solid non-linear features. The optimization algorithm used for the training of the neural network is ADAM. During training, the learning rate for the optimization algorithm is dynamically adjusted using learning rate schedulers. It is mainly decreased during the training which helps in converging towards a more accurate solution and avoids overfitting.

## 4. Model validation and discussion

To validate the performance of the hybrid TP model against the non-linear TP model, multiple validation simulations are performed. The models are solved using ODE solvers through the Julia package "DifferentialEquations.jl". A particular attention was given to code optimization to ensure type stability and to profit from code compilation. For the time-stepping algorithm, both stiff and non-stiff methods works well for the TP models. Among these, a 3-stage, 5th order Runge-Kutta method with a 3rd order embedded method for error estimation was found to be the most efficient and accurate. Therefore, this algorithm is used for all the TP models in this study.

The effect of the number of CVs used in the hybrid model on its accuracy and computational performance has been analyzed to choose the optimal number. To perform this analysis, hybrid TP model is built for different $n_{hyb}$. The performance of these hybrid models are evaluated by validating each model using the same simulation sequence used for grid convergence study of non-linear TP model (Section 2.4). Fig. 8 illustrates the change in performance of the hybrid models with different $n_{hyb}$. From this plot, it appears that reducing $n_{hyb}$ does not affect the accuracy of the model measured by the MSE until reaching a minimal number (~25 in our case), below which the accuracy deteriorates rapidly. However, as $n_{hyb}$ increases, the computational time rises significantly. In this study, the optimal number of CVs used in the linear model within the hybrid TP model is chosen as $n_{hyb} = 25$ which excels in terms of both accuracy and computational performances.

Each validation simulation consists of a number of consecutive pairs of charging and discharging processes. Similar to the training data generation process, the fluid inlet temperature is varied randomly in gradual and smooth way in all sequences, adhering to the specified ranges for charging process (0.6 to 1.0) and discharging process (0.0 to 0.2). But, for the inlet mass flow rate, even though it is varied gradually and smoothly in random aspect during the training data generation, in the validation phase it is randomly varied in a step like way. The input sequences for the validation phase are made to undergo these kind of variations to mimics the nature of these inputs in actual real world scenarios.

Each charging and discharging process lasts for a duration of $\tau = 10$. In terms of validation methodology, for each charging or discharging process, first the linear TP model is solved. Then, the corresponding linear solution is mapped to non-linear solution using the ML module. Then the final state of this non-linear solution is given as the initial condition for the next process. This is done repeatedly for all the processes in a simulation. Overall, the linear solution is mapped into
<table><caption>**Table 5**
Configuration of validation simulations.</caption><thead><tr><td>Simulation</td><td>Number of charging and discharging pairs</td><td>Total storage operation duration [h]</td><td>m Range [kg/s]</td></tr></thead><tbody><tr><td>1</td><td>4</td><td>3.3</td><td>0.45 to 0.65</td></tr><tr><td>2</td><td>10</td><td>8.3</td><td>0.45 to 0.65</td></tr><tr><td>3</td><td>20</td><td>16.6</td><td>0.45 to 0.65</td></tr><tr><td>4</td><td>30</td><td>25.0</td><td>0.45 to 0.65</td></tr><tr><td>5</td><td>30</td><td>25.0</td><td>0.35 to 0.75</td></tr><tr><td>6</td><td>40</td><td>33.3</td><td>0.35 to 0.75</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Effect of $n_{hyb}$ in the hybrid TP model's performance.</figcaption></figure>

a non linear solution at the end of each charging/discharging process. It is also possible to dynamically map the linear solution to a non-linear solution periodically, at each time step, or at a specific time. This allows to dynamically couple the hybrid TP model of PBTES with other energy system models.

The details regarding the number of pairs of processes in a simulation and the inlet mass flow rate range used can be found in Table 5. In validation simulations 1 to 4, the inlet mass flow rate is varied within the ML module's training range. The validation simulations 5 and 6 are designed to test the robustness and generality of the hybrid model, showcasing the capability of the model to perform well under different conditions which are not represented in the training data. In these simulations, the inlet mass flow rate is varied beyond the training data's scope by extending the variation bounds beyond the training range.

A detailed overview of the results of the first simulation which has only 4 pairs of charging and discharging processes is presented in the following. The dynamics of the inputs for this first validation simulation is shown in Fig. 9. For this first simulation, the two main aspects which helps in analyzing the performance of the model are plotted. Firstly, the solid temperature distributions along the length of the storage at the end of each charging process are plotted in Fig. 10. Secondly, the fluid outlet temperatures during each discharging process over respective time period are plotted in Fig. 11. In this validation study, the performance metrics such as Mean Squared Error (MSE), $R^2$ score and overall error percentage (Relative Mean Absolute Error) are used for evaluation. These metrics are calculated by comparing the temperature solutions (both fluid and solid) of the linear and hybrid TP models against the corresponding solutions from the non-linear TP model in °C. Table 6 compares the performance the metrics of the linear and hybrid TP model vs non-linear TP model, while Fig. 12 presents the MSE box plot for all the validation simulations. Table 7 provides the comparison of the computation time taken by different model to perform different validation simulations.

The hybrid TP model demonstrated significant performance in computational efficiency and accuracy. For the first validation simulation, in Figs. 10 and 11, both solid temperature distribution solution during charging processes and the fluid outlet temperature solution during the discharging processes obtained from the hybrid TP model show very

good agreement with the non-linear TP model solutions. The linear TP model, however, shows different dynamics due to the linearization process. The ML module in the hybrid model effectively maps the linear to non-linear solutions, ensuring the hybrid model's results match well with the non-linear model, demonstrating exceptional performance, as seen in these plots. From Table 6 and Fig. 12 shows that in all validation simulations, the hybrid TP model consistently achieved a lower overall MSE of 0.15 °C and an error percentage of just 0.16, and demonstrating significantly better accuracy performance than the linear TP model. An $R^2$ score of 0.9995 indicates a high level of accuracy and a strong correlation between hybrid and non-linear TP model solutions. Furthermore, the hybrid TP model also performed really good in validation simulations 5 and 6, where the inlet mass flow rate is varied beyond the training range, which also has same level of accuracy performance as other validation simulations. To more showcase these results for simulations 5 and 6, a violin plot is presented in Fig. 13, which is created by calculating the ratio of the hybrid solution to the non-linear solution at each time step and across spatial discretization. This plot shows very narrow probability distribution centered around 1 for both simulation 5 and 6. This demonstrates the hybrid model's accuracy outside the training range and highlights its explainability and generality for the scenarios outside the training range.

In terms of computational efficiency, the hybrid TP model demonstrated a remarkable improvement over the traditional non-linear TP model. In average, during the model validation as shown in Table 7, the hybrid TP model is 350 times faster than non-linear TP model. Since the linear model used in the hybrid TP model is solved using a coarse mesh with only 25 CVs, it only requires solving 50 equations (2 equations per CV). This significantly reduces computational requirements, leading to faster computation times and lower memory usage. Additionally, the ML module in the hybrid model has very low inference time, which minimally impacts the overall computational performance. Consequently, the hybrid model provides the non-linear solution much more quickly. In contrast, both the linear and non-linear TP models are typically solved for 1000 CVs, requiring 2000 equations to be solved. Specifically, the non-linear TP model is significantly slower due to the increased computational demand from its non-linearity. This requires updating more parameters for each function call and needing more iterations of solvers to maintain accuracy.

Overall, the hybrid TP model leverages the computational speed of linear models while integrating the precision of non-linear mappings through machine learning. Also the usage of fewer CVs in the hybrid TP model does not affect the accuracy and it effectively maps the coarse linear solution into high resolution non-linear solution. The hybrid TP model is highly suitable for large-scale thermal energy storage system simulations, where both rapid computation and accuracy are crucial. Since the models are studied in dimensionless form, this hybrid modeling approach ensures seamless scaling and adaptation to both small and large-scale systems. Due to its linear nature and reduced computational resource requirements, the hybrid TP model is highly compatible for optimization tasks, long-term simulations, and energy management applications. The computational efficiency of the hybrid TP model makes it an ideal choice for scenarios demanding high performance and reliability. Additionally, it allows for more simulations to be run in parallel, enabling more extensive analysis within the same time frame.
<figure><img src="image_12.png" /><figcaption>Fig. 9. Input sequences for Validation Simulation 1: (a) Inlet fluid temperature, (b) Inlet fluid mass flow rate.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Comparison of solid temperature distribution solutions between different models at the end of the 1st (a), 2nd (b), 3rd (c), and 4th (d) charging processes for the first validation simulation.</figcaption></figure>

## 5. Conclusion

In this paper a hybrid model for modeling PBTES was proposed. It combines a linearized TP model solved using a coarse FV mesh

with a ML module that maps the linear and coarse solution to a non-linear and fine solution. The proposed hybrid TP model for PBTES has demonstrated a promising gain in computational efficiency and has very good accuracy compared to the fine non-linear TP model. The hybrid TP model bridges effectively the gap between the demands of
<figure><img src="image_14.png" /></figure>

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

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

<figure><img src="image_17.png" /><figcaption>Fig. 11. Comparison of fluid outlet temperature between different models over the respective time periods during the 1st (a), 2nd (b), 3rd (c), and 4th (d) discharging processes for the first validation simulation.</figcaption></figure>

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

<figure><img src="image_19.png" /><figcaption>Fig. 12. Mean squared error bar plot for different validation sequences. (a) linear vs non-linear TP model and (b) hybrid vs non-linear TP model.</figcaption></figure>

<table><caption>**Table 6**

Comparative Accuracy Metrics: Linear and Hybrid TP vs. Non-linear TP model's solutions in °C across different validation simulations.</caption><thead><tr><th rowspan="2">Validation simulation</th><th colspan="3">Linear vs Non-linear</th><th colspan="3">Hybrid vs Non-linear</th></tr><tr><th>MSE [°C]</th><th>R<sup>2</sup></th><th>Error %</th><th>MSE [°C]</th><th>R<sup>2</sup></th><th>Error %</th></tr></thead><tbody><tr><td>1</td><td>54.19</td><td>0.6621</td><td>3.35</td><td>0.14</td><td>0.9981</td><td>0.16</td></tr><tr><td>2</td><td>43.20</td><td>0.6280</td><td>2.96</td><td>0.12</td><td>0.9989</td><td>0.14</td></tr><tr><td>3</td><td>43.17</td><td>0.6318</td><td>2.95</td><td>0.12</td><td>0.9990</td><td>0.14</td></tr><tr><td>4</td><td>44.12</td><td>0.5763</td><td>3.04</td><td>0.15</td><td>0.9992</td><td>0.16</td></tr><tr><td>5</td><td>42.63</td><td>0.6277</td><td>2.97</td><td>0.17</td><td>0.9990</td><td>0.18</td></tr><tr><td>6</td><td>43.19</td><td>0.6368</td><td>2.96</td><td>0.18</td><td>0.9990</td><td>0.18</td></tr></tbody></table>
<table><caption>**Table 7**

Computation time comparison between models for different validation sequences.</caption><thead><tr><th rowspan="2">Validation simulation</th><th colspan="3">Computation time [s]</th></tr><tr><th>Linear TP (1000 CVs)</th><th>Non-linear TP (1000 CVs)</th><th>Hybrid TP (25 CVs)</th></tr></thead><tbody><tr><td>1</td><td>0.454</td><td>20.249</td><td>0.045</td></tr><tr><td>2</td><td>1.209</td><td>40.280</td><td>0.120</td></tr><tr><td>3</td><td>2.358</td><td>83.478</td><td>0.286</td></tr><tr><td>4</td><td>3.597</td><td>125.769</td><td>0.368</td></tr><tr><td>5</td><td>3.599</td><td>127.245</td><td>0.364</td></tr><tr><td>6</td><td>4.794</td><td>182.896</td><td>0.485</td></tr></tbody></table>

<figure><img src="image_20.png" /><figcaption>**Fig. 13.** Violin plot for the ratio between hybrid and non-linear solutions for validation simulations 5 and 6.</figcaption></figure>

computational resources and the necessity for accuracy, by melding the computational speed of the coarse linear TP model with the precision of the fine non-linear TP model. In validation, the hybrid TP model is on average 350 times faster than the non-linear TP model, while also demonstrating superior accuracy, with R² score of 0.99, overall MSE of just 0.15 °C and a very low error percentage of 0.16. By also validating the model outside the training scope, model's extrapolation ability, its general applicability and robustness are shown. Thus, the hybrid TP model is capable of adapting and generating the non-linear solutions accurately under conditions across diverse scenarios that are not represented in the training phase. Furthermore, this hybrid modeling approach ensures the model's versatility, allowing it to be adapted to various storage materials and fluids, and easily scaled for both small and large-scale systems. This research opens several avenues for further exploration. One potential direction for future research could involve extending the methodologies and findings of this study to latent heat version of PBTES, which involves phase change in the filler material. Overall, the hybrid TP model is computationally less expensive and has good accuracy, which makes it highly compatible for the model optimization for control applications and simulation of PBTES over long duration which is typically the case for energy planning models.

### CRediT authorship contribution statement

**Shri Balaji Padmanabhan:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Resources, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Mohamed Tahar Mabrouk:** Writing – review & editing, Validation, Supervision, Software, Methodology, Investigation, Formal analysis, Conceptualization. **Bruno Lacarrière:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

The code used to perform the research presented in this article is available at the following link: https://gitlab.imt-atlantique.fr/hybrid_modeling/heat-storage.

### Acknowledgments

This research was conducted under the auspices of the ValaDoE chair at IMT Atlantique in partnership with Télécom Paris and Mines Saint Etienne, and was supported by Enedis, Région Pays de la Loire, Nantes Métropole and Akajoule. This work was performed using HPC resources from GENCI-IDRIS (Grant 20XX-AD011014376R1).

### References

[1] Abhishek Gautam, R.P. Saini, A review on technical, applications and economic aspect of packed bed solar thermal energy storage system, J. Energy Storage 27 (2020) 101046.

[2] Li Peng, Hongjun Wu, Qianjun Mao, Numerical and experimental study on the performance of a thermal energy storage in a solar building, J. Energy Storage 61 (2023) 106745.

[3] Edward Barbour, Dimitri Mignard, Yulong Ding, Yongliang Li, Adiabatic compressed air energy storage with packed bed thermal energy storage, Appl. Energy 155 (2015) 804–815.

[4] G. Zanganeh, A. Pedretti, A. Haselbacher, Aldo Steinfeld, Design of packed bed thermal energy storage systems for high-temperature industrial process heat, Appl. Energy 137 (2015) 812–822.

[5] Harmeet Singh, R.P. Saini, J.S. Saini, A review on packed bed solar energy storage systems, Renew. Sustain. Energy Rev. 14 (3) (2010) 1059–1069.

[6] Antoni Gil, Marc Medrano, Ingrid Martorell, Ana Lázaro, Pablo Dolado, Belén Zalba, Luisa F. Cabeza, State of the art on high temperature thermal energy storage for power generation. Part 1—Concepts, materials and modellization, Renew. Sustain. Energy Rev. 14 (1) (2010) 31–55.

[7] Hitesh Bindra, Pablo Bueno, Jeffrey F. Morris, Reuel Shinnar, Thermal analysis and exergy evaluation of packed bed thermal storage systems, Appl. Therm. Eng. 52 (2) (2013) 255–263.

[8] Thibaut Esence, Arnaud Bruch, Sophie Molina, Benoit Stutz, Jean-François Fourmigué, A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems, Sol. Energy 153 (2017) 628–654.

[9] A. Bruch, J.-Fr. Fourmigué, R. Couturier, S. Molina, Experimental and numerical investigation of stability of packed bed thermal energy storage for CSP power plant, Energy Procedia 49 (2014) 743–751.

[10] Donald E. Beasley, John A. Clark, Transient response of a packed bed for thermal energy storage, Int. J. Heat Mass Transfer 27 (9) (1984) 1659–1669.

[11] A.V. Kuznetsov, An analytical solution for heating a two-dimensional porous packed bed by a non-thermal equilibrium fluid flow, Appl. Sci. Res. 55 (1995) 83–93.

[12] Jon T. Van Lew, Peiwen Li, Cho Lik Chan, Wafaa Karaki, Jake Stephens, Analysis of heat storage and delivery of a thermocline tank having solid filler material, 2011.

[13] Tawfiq Chekifi, Moustafa Boukraa, CFD applications for sensible heat storage: A comprehensive review of numerical studies, J. Energy Storage 68 (2023) 107893.

[14] Evgeny V. Votyakov, Aristides M. Bonanos, A perturbation model for stratified thermal energy storage tanks, Int. J. Heat Mass Transfer 75 (2014) 218–223.

[15] Rocío Bayón, Esther Rojas, Analytical function describing the behaviour of a thermocline storage tank: A requirement for annual simulations of solar thermal power plants, Int. J. Heat Mass Transfer 68 (2014) 641–648.

[16] Mohamed Tahar Mabrouk, Abdelhamid Kheiri, Michel Feidt, Using generalized integral transforms to solve a perturbation model for a packed bed thermal energy storage tank, Int. J. Heat Mass Transfer 84 (2015) 633–641.

[17] Zhaoyu He, Weimin Guo, Peng Zhang, Performance prediction, optimal design and operational control of thermal energy storage using artificial intelligence methods, Renew. Sustain. Energy Rev. 156 (2022) 111977.

[18] Argyrios Anagnostopoulos, Theofilos Xenitopoulos, Yulong Ding, Panos Seferlis, An integrated machine learning and metaheuristic approach for advanced packed bed latent heat storage system design and optimization, Energy 297 (2024) 131149.
[19] Yuanji Li, Xinyu Huang, Xiaoyong Huang, Xinyu Gao, Rukun Hu, Xiaohu Yang, Ya-Ling He, Machine learning and multilayer perceptron enhanced CFD approach for improving design on latent heat storage tank, Appl. Energy 347 (2023) 121458.

[20] Kemal Ermis, Aytunc Erek, Ibrahim Dincer, Heat transfer analysis of phase change process in a finned-tube thermal energy storage system using artificial neural network, Int. J. Heat Mass Transfer 50 (15-16) (2007) 3163-3175.

[21] Z. Li, S.T. Lv, Performance analysis and optimization of packed-bed TES systems based on ensemble learning method, Energy Rep. 8 (2022) 8165-8176.

[22] Leila Darvishvand, Vahid Safari, Babak Kamkari, Meysam Alamshenas, Masoud Afrand, Machine learning-based prediction of transient latent heat thermal storage in finned enclosures using group method of data handling approach: a numerical simulation, Eng. Anal. Bound. Elem. 143 (2022) 61-77.

[23] Maarten G. Poirot, Rick H.J. Bergmans, Bart R. Thomson, Florine C. Jolink, Sarah J. Moum, Ramon G. Gonzalez, Michael H. Lev, Can Ozan Tan, Rajiv Gupta, Physics-informed deep learning for dual-energy computed tomography image processing, Sci. Rep. 9 (1) (2019) 17709.

[24] Yilin Zhuang, Yixuan Liu, Akhil Ahmed, Zhengang Zhong, Ehecatl A. del Rio Chanona, Colin P. Hale, Mehmet Mercangöz, A hybrid data-driven and mechanistic model soft sensor for estimating CO2 concentrations for a carbon capture pilot plant, Comput. Ind. 143 (2022) 103747.

[25] Ward Quaghebeur, Elena Torfs, Bernard De Baets, Ingmar Nopens, Hybrid differential equations: integrating mechanistic and data-driven techniques for modelling of water systems, Water Res. 213 (2022) 118166.

[26] Jie Li, Manu Suvarna, Lanjia Pan, Yingru Zhao, Xiaonan Wang, A hybrid data-driven and mechanistic modelling approach for hydrothermal gasification, Appl. Energy 304 (2021) 117674.

[27] Alina Kloss, Stefan Schaal, Jeannette Bohg, Combining learned and analytical models for predicting action effects from sensory data, Int. J. Robotics Res. 41 (8) (2022) 778-797.

[28] Peter Sadowski, David Fooshee, Niranjan Subrahmanya, Pierre Baldi, Synergies between quantum mechanics and machine learning in reaction prediction, J. Chem. Inf. Model. 56 (11) (2016) 2125-2128.

[29] Tianfang Xu, Albert J. Valocchi, Data-driven methods to improve baseflow prediction of a regional groundwater model, Comput. Geosci. 85 (2015) 124-136.

[30] Thibaut Esence, Arnaud Bruch, Jean-François Fourmigué, Benoit Stutz, A versatile one-dimensional numerical model for packed-bed heat storage systems, Renew. Energy 133 (2019) 190-204.

[31] Pitam Chandra, D.H. Willits, Pressure drop and heat transfer characteristics of air-rockbed thermal storage systems, Sol. Energy 27 (6) (1981) 547-553.

[32] Zhen Yang, Suresh V. Garimella, Thermal analysis of solar thermal energy storage in a molten-salt thermocline, Sol. Energy 84 (6) (2010) 974-985.

[33] Chao Xu, Zhifeng Wang, Yaling He, Xin Li, Fengwu Bai, Sensitivity analysis of the numerical study on the thermal performance of a packed-bed molten salt thermocline thermal storage system, Appl. Energy 92 (2012) 65-75.

[34] Ryan Anderson, Liana Bates, Erick Johnson, Jeffrey F. Morris, Packed bed thermal energy storage: A simplified experimentally validated model, J. Energy Storage 4 (2015) 14-23.

[35] Christopher Rackauckas, Qing Nie, DifferentialEquations.jl-a performant and feature-rich ecosystem for solving differential equations in Julia, J. Open Res. Softw. 5 (1) (2017).

[36] Uladzislau Kapustsin, Utku Kaya, Thomas Richter, A hybrid finite element/neural network solver and its application to the Poisson problem, PAMM 23 (3) (2023) e202300135.

[37] Nils Margenberg, Dirk Hartmann, Christian Lessig, Thomas Richter, A neural network multigrid solver for the Navier-Stokes equations, J. Comput. Phys. 460 (2022) 110983.

[38] Akhil Ahmed, Ehecatl Antonio del Rio-Chanona, Mehmet Mercangöz, Learning linear representations of nonlinear dynamics using deep learning, IFAC-PapersOnLine 55 (12) (2022) 162-169.

[39] Yiqi Gu, Michael K. Ng, Deep neural networks for solving large linear systems arising from high-dimensional problems, SIAM J. Sci. Comput. 45 (5) (2023) A2356-A2381.

[40] John N.R. Jeffers, Two case studies in the application of principal component analysis, J. R. Stat. Soc. Ser. C. Appl. Stat. 16 (3) (1967) 225-236.

[41] Ian Joliffe, Principal component analysis, in: Encyclopedia of Statistics in Behavioral Science, Wiley Online Library, 2005.