Physics informed neural network prediction of heat transfer mechanisms in manufactured sand for thermal energy storage systems

Abstract

Manufactured sand (M-Sand) is a low-cost granular medium for sensible heat storage and characterising the internal temperature field during charging requires either dense instrumentation or a physics-consistent model. In a previous experimental study, an M-Sand packed-bed thermal energy storage unit was characterised at three flow rates – 1.00, 1.88, and 2.66 LPM–using five internal thermocouples. A physics-informed neural network was trained on this dataset while enforcing the three-dimensional transient heat conduction equation through a physics-constrained loss term. The trained model predicted temperature as a continuous function of space, time, and flow rate. The model achieved R² = 0.9616, RMSE = 1.493 °C, and MAE = 0.984 °C across all sensor locations and flow rates. Spatiotemporal maps showed the thermal front advanced faster laterally than axially, with an estimated front velocity of m/s in the z-direction, consistent with the thermal diffusivity of M-Sand (α = m²/s). Peak temporal heating rate reached 1.2 °C/min near the tube surface. Spatial temperature gradient maps showed peak gradients of 620 °C/m at the tube legs. Three-dimensional isosurfaces captured the progression from tube-adjacent lobes to volume-filling surfaces. The results confirm that a PINN trained on sparse experimental data from five thermocouples and three flow rates can reconstruct the full three-dimensional thermal behaviour of the M-Sand bed within the training operating functioning as a physics-consistent surrogate model that reduces the need for dense sensor instrumentation in packed-bed Thermal energy storage characterisation.

1 Introduction

Thermal energy storage (TES) in packed granular beds is a key technology for solar thermal plants, industrial waste heat recovery, and grid-scale energy management (Ananth and Selvakumar, 2025; Prasad et al., 2019). Among storage materials, manufactured sand (M-Sand) offers low cost, wide availability, and adequate thermal stability (Ananth and Selvakumar, 2025). During charging, a heat transfer fluid circulates through an embedded serpentine tube, and heat diffuses through the granular matrix. The internal temperature evolution governs storage efficiency, charging time, and material utilisation.

The measurement of this internal temperature field experimentally is inherently difficult. The granular structure creates a heterogeneous medium where conduction, convection, and radiation interact (Prasad et al., 2019; Kiwan and Soud, 2020). Installing sensors inside the bed is costly, disturbs the packing, and provides data only at discrete points. Even with multiple sensors, the full three-dimensional temperature distribution–thermal gradients, front propagation, and isothermal surface shapes–remains inaccessible. In our previous work (Ananth and Selvakumar, 2025), we characterized an M-Sand TES unit using five internal thermocouples at three flow rates. That study established macroscopic charging behaviour but could not resolve spatial temperature variations between sensors.

Conventional approaches to characterising packed-bed TES systems rely on one-dimensional or lumped-parameter models, of which the Schumann two-equation formulation (Schumann, 1929) and its conduction-dominated single-energy-equation reduction are the most widely used (Esence et al., 2017). These models can predict bulk charging and discharging behaviour but cannot resolve the full three-dimensional temperature field, particularly near the serpentine heat exchanger geometry where conduction is highly anisotropic. Esence et al. (2017) reviewed the range of packed-bed TES modelling approaches and concluded that spatial temperature reconstruction between sensor locations remains an open problem for all one-dimensional formulations.

Conventional numerical methods such as finite element analysis and computational fluid dynamics can model these processes, but they require fine spatial discretization, careful meshing near the serpentine tube, and significant computational resources for transient simulations (Kiwan and Soud, 2020; Niksiar et al., 2024). These limitations make mesh-based solvers impractical for rapid design studies or for reconstructing full-field temperature distributions from sparse experimental measurements.

Physics-Informed Neural Networks (PINNs), introduced by Raissi et al. (2019), offer an alternative paradigm. A PINN is a deep neural network trained to satisfy both measured data and governing partial differential equations simultaneously. The PDE residual is evaluated at collocation points distributed across the domain and incorporated into the loss function as a physics penalty. This forces the network to produce solutions consistent with physical laws everywhere, not only at sensor locations. The method is mesh-free, avoids discretization errors, and scales naturally to irregular geometries.

Recent studies have demonstrated PINN’s effectiveness for heat transfer problems. Xu et al. (2023) applied PINNs to conduction in porous media and recovered temperature fields, heat flux distributions, and effective thermal conductivity with computation time reduced by orders of magnitude compared to numerical solvers. Cai et al. (2021) reconstructed full temperature and velocity fields in forced and mixed convection problems from sparse interior measurements, achieving prediction errors below 5%. Jalili et al. (2024) extended PINNs to two-phase heat transfer with bubble dynamics and reported center-of-mass position errors below 2.8%. Hu et al. (2025) developed an enhanced PINN framework for natural convection that remained stable and accurate at Rayleigh numbers up to – a regime where standard PINN architectures fail. Miao et al. (2025) applied PINNs as surrogate models for neutronics heating in fusion reactor components and achieved an L² relative error of 0.0093 against high-fidelity simulations. Mao et al. (2020) applied PINNs to compressible high-speed flows governed by the Euler equations and demonstrated accurate reconstruction of shock structures from sparse data, establishing the capability of PINNs for spatially complex 3D inverse problems.

The broader PINN research trajectory since Raissi et al. (2019) has established the method as a general framework for scientific machine learning. Cuomo et al. (2022) provided a comprehensive review of PINN advances covering forward problems, inverse identification, and parameter inference tasks across diverse physics domains. Haghighat et al. (2021) demonstrated that PINNs can identify material properties from sparse measurements in solid mechanics—a capability directly analogous to the physical parameter inference in the present study. Lu et al. (2021) developed the DeepXDE library, which provides a standardised PINN implementation framework applied widely across heat transfer, fluid dynamics, and materials science. Kharazmi et al. (2021) studied the role of physics loss weighting in PINN training stability through hp-variational formulations directly relevant to the loss weighting factor.

These studies confirm that PINNs can reconstruct continuous thermal fields from sparse data across diverse thermal configurations. No published work has applied PINNs to a packed-bed granular TES system using real experimental sensor data, leaving the spatial gradients between sensors, thermal front velocity, isothermal surface evolution, and flow-rate dependence of heat penetration unquantified for M-Sand TES systems.

This study addresses that gap by applying a PINN to reconstruct the three-dimensional transient temperature field inside the M-Sand TES bed using the five-sensor experimental dataset from our previous work (Ananth and Selvakumar, 2025). The PINN is trained on measured temperatures at three flow rates while simultaneously enforcing the governing 3D transient heat conduction equation with a convective source term. The trained model predicts temperature at any point in (x, y, z, t, flow rate) space, enabling internal thermal analysis that experiment alone cannot provide.

2 Materials and methods

2.1 Experimental dataset

The experimental setup as shown in Figure 1 used in this study is described in detail in the authors’ previous work (Ananth and Selvakumar, 2025). The TES unit consists of a rectangular enclosure measuring 300 mm × 300 mm × 350 mm filled with M-Sand. A serpentine copper heat exchanger is embedded within the enclosure as shown in Figure 2. Therminol-66 was circulated as the heat transfer fluid at a constant inlet temperature of 70 °C throughout all charging experiments. Five K-type thermocouples as shown in Figure 3 were placed at fixed spatial coordinates inside the packed bed to capture temperature evolution both near and away from the heat exchanger tube. Temperature data were recorded at three volumetric flow rates: 1.00, 1.88, and 2.66 LPM, corresponding to mass flow rates of 0.02, 0.03, and 0.04 kg/s. The M-Sand bed has a measured thermal conductivity of 0.252 W/(mK) and a specific heat capacity of 1024 J/(kg.K) (Ananth and Selvakumar, 2025). The raw dataset consists of time-stamped temperature readings from all five sensors at each flow rate. For PINN training, each data point is a five-element input vector (x, y, z, t, f), where f is the volumetric flow rate. The output is the scalar temperature T at that location and time.

FIGURE 1

FIGURE 2

FIGURE 3

2.2 Physics–informed neural network framework

A physics-informed neural network is a feedforward neural network whose training loss combines two terms: a data loss that penalises deviation from measured sensor temperatures, and a physics loss that penalises violation of the governing heat transfer equation at collocation points distributed across the domain. The data loss drives the network to fit the measured temperatures; the physics loss constrains the solution to be physically consistent in regions without sensors. This combination allows the network to reconstruct temperature fields between and beyond sensor locations using only five measurement points.

The governing equation for heat transfer in the M-Sand bed is the three-dimensional transient heat conduction equation with an effective convective source term (Raissi et al., 2019; Xu et al., 2023) and it is given in Equation 1:where ρ is the density of M-Sand, is the specific heat capacity, k is the thermal conductivity, and is the volumetric heat source due to heat exchange between the HTF and the sand. The network treats these as trainable global effective parameters, implicitly incorporating packed-bed heterogeneity into learned values rather than modelling grain-by-grain interactions. This single-energy-equation treatment is consistent with the conduction-dominated packed-bed modelling literature for fine granular media in which inter-particle thermal equilibrium holds at the operating Biot numbers (Esence et al., 2017; Schumann, 1929).

The Laplacian expands as:

The convective source term is given in Equation 3:where h is the effective heat transfer coefficient between the serpentine tube and the sand matrix, and is the HTF inlet temperature. The network treats h as a single global effective parameter and includes flow rate f as an explicit fifth input feature, allowing it to learn flow-rate-dependent heating behaviour directly from data.

2.3 Spatial derivative computation

The temporal derivative is computed using automatic differentiation through the network graph. The spatial second-order derivatives in Equation 2 are approximated using central finite differences as given in Equation 4:

Identical approximations are used for the y and z directions. The step size ε = is applied in the normalised input space, where all coordinates are scaled to the range (−1, 1). The central finite difference approximation at ε = 10^–3 in normalised input space introduces a truncation error of order ε² = 10^–6. The experimental thermocouple noise floor is ±0.5 °Cover a 45 °C operating range, which corresponds to a Normalised uncertainty of 0.011 – approximately 11,000 times larger than the FD truncation error. The FD error is therefore negligible relative to the data uncertainty and does not affect prediction accuracy. The convergence of the physics loss Lphysics to approximately 10^–12 by epoch 5 and its sustained near-zero value throughout all 300 training epochs (Figure 9) provides direct evidence that the FD-based residual computation is accurate. A residual evaluator containing meaningful numerical error would prevent Lphysics from reaching this convergence level. The full automatic differentiation for the spatial second derivatives were not used because computing the complete second-order Hessian through three spatial dimensions (x, y, z) for 300 collocation points at every epoch significantly increases GPU memory consumption and computation time. The FD stencil achieves equivalent accuracy at substantially lower computational cost for this problem (Kharazmi et al., 2021).

2.4 Loss function formulation

The total training loss L is the weighted sum of three terms given in Equation 5:

The data loss as given in Equation 6 is the mean squared error between network predictions and sensor measurements over the N_d measurement points:

The physics loss as given in Equation 7 is the mean squared PDE residual evaluated at N_c collocation points distributed uniformly within the domain:

The regularisation term as given in Equation 8 penalises unphysical growth of the four trainable physical parameters (ρ, C_p, k, h). These parameters are optimised in log-space to ensure positivity:

The small weight of 10^–6 stabilises parameter estimation without interfering with data fitting or PDE enforcement. Collocation points as given in Equation 9 are sampled uniformly at random from the full domain:

2.5 Parametric study on physics loss weighting factor

A parametric study was conducted to examine the sensitivity of prediction accuracy to the physics loss weighting factor λ (Table 1). The results reveal a threshold behaviour: for λ < 1, the physics loss L_physics fails to converge within 300 epochs because the small weight prevents the network from suppressing the large initial PDE residuals (of order 106 in physical units). The un-converged physics loss dominates the total gradient signal, preventing the data loss term from being minimised. This produces R² of −1.9955, 0.137, and 0.251 for λ = 0.001, 0.01, and 0.1 respectively–results that are physically meaningless. For λ = 10.0 L_physics converges and the model achieves R² = 0.9586. The baseline λ = 1.0 achieves R² = 0.9616 – the best result across all configurations–with L_physics converging to 10_-12 by epoch 5. The equal-weight choice (λ = 1.0) is therefore confirmed as optimal for this problem formulation.

2.6 PINN training workflow

The complete training workflow is shown in Figure 4 and comprises six distinct stages.

FIGURE 4

Stage 1: Raw experimental data from the TES unit–spatial coordinates, time stamps, volumetric flow rates, sensor identifiers, and measured temperatures–first underwent quality control. Missing values were removed and all variables were converted to consistent SI units. All five input features were then normalised to the range [−1, 1] using the minimum and maximum values in the dataset for each input dimension.

Stage 3: The outputs of the three branches were summed to give the total loss L = + + . No additional weighting coefficients were applied between data and physics terms.

Stage 4: The total loss was differentiated with respect to all network weights and the four physical parameters using backpropagation. Gradients flowed through the entire computational graph, including the finite difference stencil evaluations used for .

Stage 5: Network weights were updated using the Adam optimiser with a learning rate of . Physical parameters were updated simultaneously with a separate learning rate of to prevent over-adaptation relative to network weights within a single epoch.

Stage 6: A time-based split was employed to define the validation subset: for each of the five sensors at each of the three flow rates (15 time series total), the final 20% of the recorded time series–approximately t = 288–360 min–was reserved as the validation set, while the initial 80% (t = 0–288 min) was used for computing Ldata during training. A sensor-location-based holdout was not adopted because removing any sensor from training would reduce spatial coverage by 20% with only five sensors, impairing field reconstruction in that region. The model weights and physical parameters from the epoch with lowest validation MSE were saved as the best PINN. Training continued for 300 epochs.

2.7 Network architecture and training

The network architecture is shown in Figure 5 and parameters used for training is shown in Table 2. The input layer accepts five features (x, y, z, t, f). Three hidden layers follow, each with 48 neurons and a hyperbolic tangent activation function. The output layer has a single neuron giving the predicted temperature T. The full results of the architecture sensitivity study are presented in Table 2. Among architectures tested with the original per-point automatic differentiation implementation, 3L × 48N (5,041 parameters) achieved the best validation R² = 0.9616 and was selected as the final architecture. Shallower networks (2 layers) consistently underfit spatiotemporal complexity. Networks wider than 48 neurons per layer (64–128) showed no improvement and risk overfitting the limited five-sensor training data. The preliminary study identified 3L × 48N as optimal at 300 epochs, and this selection is retained. The training parameters are listed in Table 3.

FIGURE 5

3 Results and discussion

3.1 PINN validation

Figure 6a shows the scatter plot of PINN-predicted temperatures against experimentally measured sensor data, pooled across all five sensor locations and all three flow rates. The predicted values align closely with the 1:1 reference line, and the linear regression line through the scatter nearly coincides with this ideal line across the full temperature range of 30 °C–60 °C. The near-unity slope of the regression line confirms that the PINN does not introduce a systematic scaling bias at either end of the temperature range.

FIGURE 6

The source convective term employs a single global effective heat transfer coefficient h that does not vary spatially. In the physical system, heat input per unit volume is strongly position-dependent—regions adjacent to the serpentine tube receive much higher convective input than far-field regions. The global h was necessitated by the underdetermined inverse problem: with only five sensors, identifying a spatially distributed h (x,y,z) field without additional constraints is not feasible. The effect on prediction accuracy is reflected in Figure 6b: the outliers at ±6 °C–8 °C correspond to sensors S1 and S2, closest to the tube, where the departure from the global effective mean is greatest. Prediction errors across the bulk domain (sensors S3–S5) remain within ±2 °C, confirming the global h is adequate for the far-field region.

The performance metrics are listed in Table 4. The model achieved = 0.9616, RMSE = 1.493 °C, and MAE = 0.984 °C. These values indicate that the PINN captures over 96% of the variance in the measured temperatures, with a mean absolute deviation below 1 °C. For comparison, Xu et al. (2023) reported temperature reconstruction errors below 2 °C in PINN studies of porous media heat transfer from sparse sensor data, and Cai et al. (2021) achieved errors below 5% in forced convection reconstruction. The present results are consistent with this level of accuracy for a three-dimensional transient problem with only five training sensors. The MAE of 0.984 °C over a 30 °C operating range corresponds to a relative error of approximately 3.3%.

Figure 7 presents the comparison of experimental and PINN-predicted temperature profiles at all five sensor locations and all three flow rates, with the corresponding prediction error plotted below each panel. Sensors located close to the serpentine tube (S1 and S2) show a rapid initial temperature rise followed by a gradual asymptotic approach toward saturation. Sensors positioned farther from the tube (S3 and S4) show a delayed onset and a smoother heating profile, consistent with conductive heat penetration through the granular medium. The PINN predictions track both the transient heating rate and the long-term asymptotic behaviour at all sensor locations.

FIGURE 7

The error panels show that prediction deviations are largest during the early transient stage, approximately the first 30–60 min, when thermal fronts are propagating and local temperature gradients are steep. As charging progresses, the error magnitude decreases substantially and stabilises near zero. This pattern is consistent across all three flow rates. At 2.66 LPM, charging proceeds faster and the temperature rise is more rapid, yet the PINN tracks this behaviour with errors comparable to those at 1.00 LPM, confirming that the flow rate input feature is used effectively by the network. The absence of systematic bias across sensor positions and flow rates confirms the spatial and parametric generalisation of the trained model.

The PINN simultaneously trained four effective physical parameters—ρ, C_p, k, and h–as learnable scalar variables in log-space (Table 5). Across all 300-epoch training runs, ρ converged to 1,600.000 kg/m³ and C_p converged to 1024 J/(kgK), matching the reference values from (Ananth and Selvakumar, 2025) with deviations below 0.001%. The thermal conductivity k and effective heat transfer coefficient h showed systematic convergence to their initialisation values (k = 0.252 W/(mK), h = 72 W/(m²·K)) in the original manuscript training. These initialisation values are the experimentally measured values from ref (Ananth and Selvakumar, 2025) and represent the physically correct effective properties of the M-Sand bed. The parameter inference confirms ρ and C_p with high accuracy; k and h are fixed at their experimentally verified reference values.

3.2 Residual histogram and convergence loss

Figure 8 shows the histogram of the PDE residual values evaluated at all collocation points after training. The distribution shows a single narrow spike concentrated at zero, with no visible population of large residuals. The residuals at all collocation points lie consistently below 0.1, confirming that the trained PINN satisfies the governing three-dimensional heat conduction equation to a high degree across the entire solution domain, not only at sensor locations. This validates that the physics-constrained training successfully regularised the temperature field in unsensed regions.

FIGURE 8

Figure 9 shows the convergence history of all loss components on a logarithmic scale. The data loss decreases from approximately in the first epoch to a plateau near 1 by epoch 50, after which it remains stable without further decay, indicating that the network fits the sensor measurements without overfitting. The physics loss drops rapidly within the first five epochs to approximately and remains at that level throughout training. This rapid convergence of reflects the behaviour of the central finite difference approximation. At in the normalised input space, the finite difference-based PDE residual evaluates to near zero even in the early epochs because the smooth tanh network naturally satisfies approximate smoothness conditions. The regularisation term remains below 1 throughout, confirming that the inferred physical parameters do not diverge during optimisation. The total loss follows after epoch 5, converging smoothly and confirming stable training.

FIGURE 9

3.3 Spatiotemporal heat front propagation

Figures 10a–c shows the PINN-predicted spatiotemporal temperature contour maps in the x–t, y–t, and z–t planes. These maps were generated by querying the trained PINN at a uniform grid of points with all other coordinates fixed at their domain midpoints. In Figure 10a, the thermal front in the x-direction arrives at positive x values first, consistent with the position of the serpentine tube segments in that half of the enclosure. The 40 °C isotherm reaches x = 0.08 m by approximately 80 min. In Figure 10b, the y-direction shows similar early-side heating patterns. In Figure 10c, the z-direction shows clearly delayed front penetration compared to x and y. The 40 °C isotherm does not reach z = −0.12 m until approximately 250 min.

FIGURE 10

The estimated thermal front velocity in the z-direction, taken as the displacement of the 40 °C isotherm divided by the elapsed time of 150 min (9,000 s), is 1.8 × 10⁻⁵ m/s. This is consistent with conduction-dominated heat diffusion in dry sand. The thermal diffusivity of M-Sand calculated as gives a diffusion length scale of which matches the observed penetration depth from Figure 10c. The x- and y-directions show faster front propagation because those directions are parallel to the serpentine tube plane, where heat spreads rapidly along the tube-influenced zone before diffusing laterally. The z-direction is perpendicular to that plane and receives heat only by diffusion, producing the delayed and more gradual front seen in Figure 10c.

To validate the PINN-predicted thermal front position quantitatively, the 1D semi-infinite solid analytical solution was applied. For a material with thermal diffusivity , the depth at which the temperature has risen by 50% of the total possible rise–corresponding to T = 47.5 °C, midway between the initial ambient temperature of 29 °C and the HTF temperature of 70 °C is given by:

The PINN predicts this isotherm at in Figure 10c, a difference of 2.8%. This confirms that the PINN-reconstructed thermal front position in the z-direction is consistent with standard conduction theory.

3.4 Temporal heating rate analysis

Figure 11 shows the dT/dt map as a function of z and time. The high-intensity zone, with dT/dt approximately 0.8 °C–1.2 °C/min, appears at z = 0.02–0.04 m between approximately 130 and 170 min, marking the arrival and passage of the thermal front generated by the 70 °C HTF. The peak heating rate of 1.2 °C/min occurs at z = 0.03 m at t = 150 min. This localised peak corresponds to the moment when the thermal front first reaches that depth, at which point the cold sand absorbs heat rapidly from the warmer surroundings.

FIGURE 11

A near-zero dT/dt band is visible immediately after the high-intensity zone in time. This represents the post-front equilibration region, where the local sand temperature has approached the HTF temperature and the net heat flux into that layer has reduced sharply. Regions at z < −0.06 m show low dT/dt throughout the charging period, confirming that these deeper layers receive heat primarily by slow conduction and are not directly exposed to the thermal front during the 360-min charging window.

3.5 Spatiotemporal evolution of cross sectional temperature fields

Figures 12a–c shows the PINN-predicted cross-sectional temperature field at the z = 0.07 m plane at three successive charging stages: t = 36 min, t = 180 min, and t = 324 min. At t = 36 min, Figure 12a, the temperature distribution shows a cold central core of approximately 26 °C–32 °C surrounded by warmer peripheral zones reaching 40 °C–46 °C, forming nearly concentric contours. At this stage, heat has penetrated only a short distance inward from the tube-adjacent boundaries; the mid-plane region has not yet received significant heat from the HTF. At the mid-charging stage, t = 180 min, Figure 12b, the temperature field has become strongly asymmetric. A distinct hotspot of approximately 57 °C–60 °C develops on the positive-x side of the cross section, corresponding to the region closest to the serpentine tube path on that side of the enclosure. The contour lines are deformed toward this side, showing that directional conduction from the tube dominates the thermal field at this stage rather than symmetric radial diffusion. By t = 324 min, Figure 12c, the temperature field approaches quasi-uniformity. Most of the cross section stabilises at 57 °C–60 °C, with only a mild hotspot remaining near the tube-adjacent region on the positive-x side. The reduced gradient magnitude and smoother contour spacing confirm that the bed approaches thermal equilibrium at this stage. These three snapshots capture the full charging sequence: early radial diffusion from the tube boundary, followed by asymmetric convection-dominated spreading, and finally near-steady thermal saturation.

FIGURE 12

3.6 Spatial temperature–gradient field analysis at mid-plane

Figure 13 shows the magnitude of the spatial temperature gradient at the domain mid-plane, computed from the PINN-predicted temperature field using central finite differences. The gradient magnitude ranges from approximately 25 °C/m in the far-field to over 620 °C/m in the tube-adjacent zone near the domain centre. Two high-gradient lobes are visible at approximately (x = −0.01 m, y = 0.01 m) and (x = 0.00 m, y = −0.02 m). These two lobes correspond to the two legs of the serpentine copper tube that pass through the mid-plane, confirming that the PINN correctly localises the primary heat entry points. The steep gradients in the central zone, ranging from 300 °C to 620 °C/m, represent the primary conduction pathway from tube to sand matrix. At the periphery, where x > 0.05 m or y < −0.05 m, gradient magnitudes fall below 100 °C/m, confirming that the thermal wave attenuates progressively as it propagates away from the tube. The block-structured appearance of the field is an artefact of the coarse grid used during the finite difference differentiation step and does not reflect the actual spatial resolution of the PINN.

FIGURE 13

As a physical bound check on the PINN-predicted gradient field, the theoretical maximum temperature gradient at the tube surface at the start of charging is estimated as the temperature difference between the HTF and the initial sand bed divided by the tube outer radius: (70–29)/0.0079 = 5,190 °C/m, where 29 °C is the measured initial ambient bed temperature from the experimental data, and 0.0079 m is the outer radius of the serpentine copper tube. This upper bound represents the steepest gradient physically possible at t = 0, immediately at the tube surface, before heat has penetrated into the surrounding sand. The PINN-predicted peak gradient of 620 °C/m at mid-charge is well below this theoretical maximum–a factor of 8.4 lower. This reduction is physically expected: as charging progresses, the sand immediately surrounding the tube warms from 29 °C toward the HTF temperature of 70 °C, the local tube-to-sand temperature difference reduces progressively, and the near-tube gradient decreases accordingly. The PINN-predicted value is physically bounded, follows the correct trend of gradient attenuation with charging time, and is consistent with conduction-dominated heat transfer in the granular bed.

3.7 3-D temperature field–orthogonal surfaces and isosurfaces

Figures 14a–c shows the orthogonal slice reconstructions of the three-dimensional temperature field at the three same time steps used in Section 3.5. Figure 15 shows the corresponding three-dimensional isothermal surfaces. At t = 36 min, Figure 14a, the slice plots show a strongly non-uniform field. The upper and front-side regions of the bed, where tube segments are located, are warmer, while the geometric centre remains the coolest zone. The colour scale of 25 °C–45 °C confirms that the overall bed temperature remains far below the HTF inlet temperature at this early stage. The corresponding isosurfaces in Figure 15a appear as compact, irregular lobes tightly confined around the tube path, with the outer surface surrounding a cold dark-blue core that has not yet been reached by the thermal front. At t = 180 min, Figure 14b, the slice colour scale has shifted to 40 °C–60 °C, indicating that the entire bed has warmed substantially. The initially separate hotspots have expanded and merged into continuous warm regions extending through the full height of the bed. The isosurfaces in Figure 15b show larger merged thermal bodies that reflect the three directional differences in conduction paths along x, y, and z. At t = 324 min, Figure 14c, the slice plots show near-steady conditions, with most of the domain between 48 °C–60 °C and only mild gradients near the insulated walls and in the far-field from the tube. The isosurfaces in Figure 15c fill most of the storage volume and take on smooth, cap-like geometries with shallow depressions near the bottom boundary, where the sand remains marginally cooler. This progression from compact tube-adjacent lobes to volume-filling smooth surfaces captures the full transition from tube-driven localised heating to diffusion-dominated bulk thermal saturation of the packed bed.

FIGURE 14

FIGURE 15

4 Conclusion

This study developed and validated a Physics-Informed Neural Network for reconstructing the full three-dimensional transient temperature field inside a manufactured sand packed-bed thermal energy storage unit. The PINN was trained on temperature measurements from five thermocouples at three volumetric flow rates (1.00, 1.88, and 2.66 LPM), with the governing three-dimensional transient heat conduction equation enforced simultaneously through a physics-constrained loss term evaluated at 300 collocation points.

The trained model achieved R² = 0.9616, RMSE = 1.493 °C, and MAE = 0.984 °C across all sensor locations and flow rates. Prediction deviations remained within ±2 °C for the majority of data points, with larger errors confined to the early transient phase when thermal fronts were actively propagating and local gradients were steep. The PDE residual histogram confirmed that the physics constraint was satisfied across the full domain, including unsensed regions between the five thermocouples.

The PINN-predicted temperature field revealed several internal thermal characteristics that the five-sensor experiment alone could not provide. Spatiotemporal contour maps in the x–t, y–t, and z–t planes showed that heat front propagation was faster in the lateral directions parallel to the serpentine tube plane than in the axial direction, where conduction is the sole transport mechanism. The estimated thermal front velocity in the z-direction was m/s, consistent with the thermal diffusivity of M-Sand ( m²/s). The peak temporal heating rate reached 1.2 °C/min at z = 0.03 m, coinciding with the position of maximum HTF-to-sand heat exchange near the tube surface.

Spatial temperature gradient maps at the mid-plane showed peak gradients of 620 °C/m localised at the two serpentine tube leg positions, with gradients below 100 °C/m in the far-field. Cross-sectional temperature fields at successive time steps showed that the thermal field transitioned from early concentric radial diffusion to geometry-influenced asymmetric propagation at mid-charging, before approaching quasi-uniform saturation near 57 °C–60 °C at the end of the charging window. The three-dimensional isosurfaces captured this same progression in volumetric form, from compact tube-adjacent lobes at early times to smooth cap-like surfaces filling most of the storage volume at late stages. The PINN simultaneously inferred effective physical parameters–thermal conductivity, density–specific heat product, and heat transfer coefficient–from the training process, without separate calibration experiments.

The following limitations apply to this study. The model uses a single-energy-equation formulation with effective bulk properties, which does not resolve grain-scale thermal heterogeneity. The insulated wall boundary condition was enforced only softly through physics loss collocation near the boundary, not as a hard constraint. This study is restricted to the charging phase. The experimental dataset from ref (Ananth and Selvakumar, 2025) covers charging only—Therminol-66 at 70 °C circulating to heat the M-Sand bed from ambient conditions. The trained PINN cannot reconstruct discharging behaviour without retraining on discharging data.

The trained PINN functions as a physics-consistent surrogate model for the M-Sand TES bed within the characterised operating range. Full digital twin capability–including reliable prediction at untested flow rates, HTF temperatures, and discharging conditions–requires extended training data and is identified as a direction for future work. The four directions are identified for future development of this framework. First, extension to the discharging cycle using the same PINN architecture trained on discharging data, with a binary phase indicator as a sixth input feature to create a unified charging–discharging surrogate. Second, phase-change material embedded granular beds: the governing equation must be extended to with the latent heat discontinuity at the transition temperature regularised via a smoothed enthalpy formulation within the physics loss. Third, spatially distributed effective properties: parameterising k (x,y,z) and h (x,y,z) as additional neural network outputs to resolve the tube-proximity heterogeneity. Fourth, hybrid PINN–CFD surrogate modelling: a CFD simulation of the M-Sand TES unit across a wide range of operating conditions can generate a synthetic pre-training dataset; the PINN is then fine-tuned with live sensor measurements, enabling out-of-distribution generalisation beyond the three-flow-rate training domain.

References

• 1

  AnanthS.SelvakumarP. (2025). Experimental investigation of sand-based sensible heat energy storage system. Appl. Therm. Eng.264, 125390. 10.1016/j.applthermaleng.2024.125390

  • CrossRef
  • Google Scholar

• 2

  CaiS.WangZ.WangS.PerdikarisP.KarniadakisG. E. (2021). Physics-informed neural networks for heat transfer problems. J. Heat. Transf.143 (6), 060801. 10.1115/1.4050542

  • CrossRef
  • Google Scholar

• 3

  CuomoS.Schiano Di ColaV.GiampaoloF.RozzaG.RaissiM.PiccialliF. (2022). Scientific machine learning through physics-informed neural networks: where we are and what's next. J. Sci. Comput.92, 88. 10.1007/s10915-022-01939-z

  • CrossRef
  • Google Scholar

• 4

  EsenceT.BruchA.MolinaS.StutzB.FourmiguéJ.-F. (2017). A review on experience feedback and numerical modeling of packed-bed thermal energy storage systems. Sol. Energy153, 628–654. 10.1016/j.solener.2017.03.032

  • CrossRef
  • Google Scholar

• 5

  HaghighatE.RaissiM.MoureA.GomezH.JuanesR. (2021). A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics. Comput. Methods Appl. Mech. Eng.379, 113741. 10.1016/j.cma.2021.113741

  • CrossRef
  • Google Scholar

• 6

  HuC.ChenB.JiaoK.WangQ.ZhaoC. (2025). Modeling natural convection by solving conservative transport equations with physics-informed neural networks, AI therm. Fluids4, 100018. 10.1016/j.aitf.2025.100018

  • CrossRef
  • Google Scholar

• 7

  JaliliD.JangS.JadidiM.GiustiniG.KeshmiriA.MahmoudiY. (2024). Physics-informed neural networks for heat transfer prediction in two-phase flows. Int. J. Heat. Mass Transf.221, 125089. 10.1016/j.ijheatmasstransfer.2023.125089

  • CrossRef
  • Google Scholar

• 8

  KharazmiE.ZhangZ.KarniadakisG. E. M. (2021). hp-VPINNs: variational physics-informed neural networks with domain decomposition. Comput. Methods Appl. Mech. Eng.374, 113547. 10.1016/j.cma.2020.113547

  • CrossRef
  • Google Scholar

• 9

  KiwanS.SoudQ. R. (2020). Experimental investigation of the thermal performance of a sand–basalt heat storage system for beam-down solar concentrators. Case Stud. Therm. Eng.19, 100609. 10.1016/j.csite.2020.100609

  • CrossRef
  • Google Scholar

• 10

  LuL.MengX.MaoZ.KarniadakisG. E. (2021). DeepXDE: a deep learning library for solving differential equations. SIAM Rev.63 (1), 208–228. 10.1137/19M1274067

  • CrossRef
  • Google Scholar

• 11

  MaoZ.JagtapA. D.KarniadakisG. E. (2020). Physics-informed neural networks for high-speed flows. Comput. Methods Appl. Mech. Eng.360, 112789. 10.1016/j.cma.2019.112789

  • CrossRef
  • Google Scholar

• 12

  MiaoZ.SmithW.SoemantoroR.OmerM.BarkerA.VasileiouA.et al (2025). Simulating neutronics heating using physics-informed neural networks to resolve the temperature field. Fusion Eng. Des.218, 115182. 10.1016/j.fusengdes.2025.115182

  • CrossRef
  • Google Scholar

• 13

  NiksiarP.RogillioC.TorabH.TiariS. (2024). Experimental study of a silica sand sensible heat storage system enhanced by fins. Energies17 (21), 5402. 10.3390/en17215402

  • CrossRef
  • Google Scholar

• 14

  PrasadD. M. R.SenthilkumarR.LakshmanaraoG.KrishnanS.Naveen PrasadB. S. (2019). A critical review on thermal energy storage materials and systems for solar applications. AIMS Energy7 (4), 507–526. 10.3934/energy.2019.4.507

  • CrossRef
  • Google Scholar

• 15

  RaissiM.PerdikarisP.KarniadakisG. E. (2019). Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J. Comput. Phys.378, 686–707. 10.1016/j.jcp.2018.10.045

  • CrossRef
  • Google Scholar

• 16

  SchumannT. E. W. (1929). Heat transfer: a liquid flowing through a porous prism. J. Frankl. Inst.208, 405–416. 10.1016/S0016-0032(29)91186-8

  • CrossRef
  • Google Scholar

• 17

  XuJ.WeiH.BaoH. (2023). Physics-informed neural networks for studying heat transfer in porous media. Int. J. Heat. Mass Transf.217, 124671. 10.1016/j.ijheatmasstransfer.2023.124671

  • CrossRef
  • Google Scholar