A two-dimensional rectangular enclosure filled with a phase change material (PCM) is considered. The enclosure has a height H and width W, with a finite reference depth D introduced for energy scaling. The PCM is initially at a uniform temperature T0, which is lower than the melting temperature Tm. n-Octadecane was selected because it is a well-documented paraffin PCM with a melting temperature close to low-temperature thermal storage and building-energy applications. Its thermophysical properties are widely available, and benchmark numerical and experimental data exist for natural-convection melting in rectangular enclosures. This makes it suitable for isolating the effect of renewable intermittency without introducing uncertainty from poorly characterized material behavior. At time t=0, the left wall is subjected to a time-dependent heat flux representing renewable energy input. The right wall is exposed to convective heat extraction, while the remaining boundaries are assumed to be adiabatic. This configuration enables direct coupling between phase-change dynamics and time-dependent energy delivery.

The formulation is based on the following physical assumptions. The flow is laminar and incompressible, and the liquid PCM is treated as a Newtonian fluid. Thermophysical properties are assumed constant, except for density variations in the buoyancy term, which are modeled using the Boussinesq approximation. Radiative heat transfer is neglected. Phase change is assumed to occur within a finite but narrow temperature interval.

Figure 1 illustrates the physical configuration of the enclosure and the imposed boundary conditions. The left boundary represents intermittent renewable heat input, the right boundary represents energy extraction, and the remaining boundaries are thermally insulated. This configuration forms the basis for the governing equations described in the following section.

The phase-change process is modeled using the enthalpy–porosity formulation, which allows simultaneous resolution of fluid flow and heat transfer on a fixed grid.

Figure 2 illustrates the overall structure of the proposed framework. The formulation integrates stochastic boundary conditions, physics-based phase-change modeling, energy-based performance evaluation, and uncertainty quantification within a unified workflow.

As shown in Figure 2, the framework begins with the definition of time-dependent boundary conditions, including stochastic renewable input and prescribed demand profiles. These inputs are processed through a physics-resolved enthalpy–porosity model, which yields the transient evolution of thermal and flow fields. The resulting fields are then transformed into energy-based quantities, enabling direct comparison between available energy and demand. Finally, reliability metrics and statistical analysis are applied to quantify system performance across multiple realizations.

The continuity equation is given by:

The momentum equation is expressed as:

The energy equation is written in terms of enthalpy:

The total enthalpy is defined as:

where the sensible enthalpy is:

and the latent enthalpy is:

The liquid fraction is defined as:

The damping term used to suppress velocity in the solid and mushy regions is given by:

where is introduced to prevent division by zero.

The renewable heat input is imposed as a time-dependent heat flux:

where is a stochastic process representing temporal fluctuations and σ is the fluctuation intensity.

Convective heat extraction at the right boundary is modeled as:

The no-slip condition is imposed on all solid walls:

Figure 3 shows the temporal characteristics of the imposed boundary conditions. The renewable input exhibits a diurnal pattern modulated by stochastic fluctuations, while the demand profile is phase-shifted relative to the input. This deliberate misalignment introduces periods of energy surplus and deficit, which are essential for evaluating system reliability.

The characteristic dimensionless numbers governing the system behavior are defined as follows:

Rayleigh number:

Stefan number:

Prandtl number:

These parameters characterize the relative importance of buoyancy-driven flow, latent heat effects, and momentum diffusivity.

The governing equations are solved using the finite-volume method on a structured and uniformly spaced grid. A collocated arrangement is adopted to ensure consistent coupling between velocity, pressure, and temperature fields.

Convective terms are discretized using a second-order upwind scheme, providing a balance between numerical stability and accuracy. Diffusive terms are discretized using second-order central differencing. Temporal discretization is performed using a first-order implicit scheme, which ensures unconditional stability for the transient phase-change problem.

The phase-change process is modeled using the enthalpy–porosity approach, in which the mushy region is treated as a porous medium governed by the permeability parameter Cmush, as listed in Table 1.

Convergence is monitored using normalized residuals for continuity, momentum, and energy equations, with thresholds of 10⁻³, 10⁻⁸, and 10⁻¹⁵, respectively. In addition, a global energy-based convergence criterion is imposed by ensuring that the variation in total stored energy between successive iterations remains negligible.

To ensure numerical reliability, grid-independence and time-step sensitivity analyses are performed. The selected discretization parameters represent a compromise between computational efficiency and accuracy, with deviations relative to refined configurations remaining within acceptable limits.

The numerical simulations were performed using an in-house finite-volume code developed specifically for phase-change heat transfer and natural-convection melting problems. The solver was implemented in MATLAB and employs a pressure-based formulation on a structured collocated grid. The use of computational fluid dynamics has become increasingly important for resolving coupled heat-transfer processes in advanced thermal-energy systems, particularly where transient flow, heat recovery, and system-level performance interact (Riffat and Samaei 2026). Pressure–velocity coupling was achieved using the SIMPLE algorithm, while convective and diffusive fluxes were discretized using second-order spatial schemes. Temporal integration was performed using an implicit first-order scheme to ensure numerical stability during the transient phase-change process.

The computational framework incorporates the enthalpy–porosity method to model the liquid–solid transition within the PCM domain and includes dedicated routines for energy tracking, reliability assessment, and stochastic renewable-input generation. All simulations were executed under identical numerical settings, and the complete workflow, including boundary-condition generation, solution procedures, post-processing, and statistical analysis, was automated to ensure consistency and reproducibility across all realizations. Pressure–velocity coupling is achieved using the SIMPLE algorithm.

Spatial discretization follows the schemes described in Section 2.5. Under-relaxation factors are applied to enhance numerical stability and ensure convergence of the coupled equations.

At each time step, convergence is achieved when the residual criteria are satisfied and the relative change in total stored energy falls below 10⁻⁵. The adequacy of the numerical setup is verified through validation, grid-independence, time-step sensitivity, and energy conservation analyses presented in Section 4.

All simulations are conducted using a version-controlled computational workflow based on the finite-volume solver described above. The reference configuration, including geometry, material properties, boundary conditions, stochastic forcing parameters, and numerical settings, is fully specified in Table 1.

A fixed random-seed registry is used for stochastic realizations to ensure repeatability. The numerical case files, input forcing profiles, output datasets, and post-processing scripts are available from the corresponding author upon reasonable request. The value of 103.13 kJ should not be interpreted as a full-scale building-storage capacity. It corresponds to the selected laboratory-scale/reference enclosure volume used for physics-resolved numerical testing. The objective was to isolate the effect of intermittency under controlled equal-energy forcing, rather than to size a practical full-scale TES unit. Because the governing response is evaluated through normalized reliability metrics, the conclusions concern sensitivity to temporal variability rather than absolute storage capacity.

The proposed framework is developed to evaluate PCM-based thermal energy storage systems from a system-level perspective, in which the primary objective is not only energy storage but reliable energy delivery under time-dependent operating conditions.

The methodology consists of two tightly coupled components. The first component is a physics-based simulation layer that resolves the phase-change process using the enthalpy–porosity method. The second component is an energy-based evaluation layer that quantifies the ability of the system to meet time-varying demand. The coupling between these components enables direct translation of internal thermal dynamics into system-level performance metrics.

Each simulation produces time-dependent fields of temperature, velocity, and liquid fraction. These fields are post-processed to obtain energy-related quantities that directly reflect the interaction between energy supply, storage, and demand. The framework operates over a full daily cycle, allowing the temporal interaction between intermittent input and time-shifted demand to be explicitly captured.

To move beyond temperature-based evaluation, system performance is formulated in terms of energy flows and demand satisfaction.

The instantaneous output power is defined as:

The cumulative available energy is obtained by time integration:

In this formulation, available energy is defined as the cumulative extractable energy delivered through the extraction boundary under the imposed heat-transfer condition. It does not represent the total internal energy stored in the PCM, but rather the portion that can be effectively supplied to the load.

The time-dependent demand profile is defined as:

The total demand over one cycle is given by:

To quantify the system’s inability to meet demand, the concept of energy deficit is introduced.

The energy deficit is defined as:

The supplied energy is defined as:

This formulation ensures that only the energy that can actually be delivered to the load is counted, thereby providing a physically meaningful measure of system performance.

Based on the energy quantities defined above, system-level reliability metrics are introduced.

The Energy Reliability Index (ERI) is defined as:

The Load Matching Index (LMI) is defined as:

Because the energy-based failure probability is directly complementary to the Energy Reliability Index, it was not retained as a primary metric in the revised framework to avoid redundant interpretation. The analysis therefore focuses on four non-redundant indicators: the Energy Reliability Index, the Load Matching Index, the cumulative energy deficit, and the critical duration of unmet demand.

The critical duration of unmet demand is defined as:

where is an indicator function equal to one when the condition is satisfied and zero otherwise.

The renewable heat input is modeled as a deterministic–stochastic hybrid signal to represent both diurnal variation and short-term fluctuations.

The imposed heat flux is defined as:

The stochastic component ξ(t)is generated using an Ornstein–Uhlenbeck process:

where is the correlation time, is the noise intensity, and is the increment of a Wiener process.

The parameter controls the macroscopic level of intermittency imposed at the boundary, whereas governs the internal generation of the stochastic signal. In the present implementation, is calibrated such that the resulting fluctuations match the prescribed intermittency level.

Three representative forcing scenarios are considered to isolate the effect of temporal variability. Scenario I corresponds to deterministic input without fluctuations, Scenario II represents moderate intermittency, and Scenario III represents strong intermittency. The total input energy over the full cycle is maintained constant across all scenarios to ensure that the observed differences arise solely from temporal variability.

Figure 4 illustrates the interaction between energy input and demand, showing that temporal misalignment between input and demand results in alternating periods of energy surplus and deficit, even when the total daily input energy remains unchanged.

The evaluation is performed over a complete daily cycle, ensuring that both charging and discharging processes are captured.

The procedure consists of solving the coupled phase-change problem, extracting the instantaneous output power, integrating the available energy, comparing it with the demand profile, and computing the resulting deficit and reliability metrics. This sequential structure ensures full traceability between physical processes and performance indicators.

To ensure statistical robustness, the stochastic simulations are repeated over an ensemble of realizations. The ensemble size is selected based on convergence of the running mean and standard deviation of the key metrics.

Table 2 presents the statistical convergence of the ensemble results, showing the evolution of ERI, LMI, and energy deficit as the number of realizations increases.

Convergence is considered achieved when the relative change in the running mean of ERI remains below 1% over the final portion of the ensemble.

To isolate the effect of intermittency and quantify the contribution of the proposed framework, a set of baseline cases is defined.

Table 3 summarizes the baseline configurations used in the study.

Case B1 represents deterministic heating without demand coupling, Case B2 introduces demand under deterministic conditions, and Case B3 corresponds to the full intermittent-forcing configuration with reliability evaluation. This structured comparison allows clear identification of the additional insight provided by the reliability-oriented framework.

The numerical formulation was validated against the benchmark melting problem for n-octadecane in a rectangular enclosure reported by Vogel and Thess (2019). The selected benchmark represents natural-convection-dominated melting and provides reference values for the global liquid fraction, mean temperature, and total heat-transfer rate at different time instances. This comparison was used to assess whether the present enthalpy–porosity model can reproduce both the early conduction-dominated stage and the later convection-dominated melting regime.

The global liquid fraction was calculated as:

The total heat-transfer rate was evaluated as:

Table 4 compares the present numerical results with the benchmark solution. The comparison includes global liquid fraction, mean temperature, and total heat-transfer rate at 2 h, 4 h, 8 h, and 12 h.

The validation results show that the relative error remains below 4% for all evaluated quantities. The liquid fraction and mean temperature follow the benchmark trends closely, while the heat-transfer rate remains consistent with the reference values throughout the melting process. These results indicate that the model captures the dominant thermal and convective mechanisms governing n-octadecane melting.

Figure 5 shows the validation of the present enthalpy–porosity model against the benchmark solution, including liquid fraction evolution, thermal response, heat-transfer behavior, and quantitative comparison of key variables.

As shown in Figure 5, the numerical model reproduces the benchmark response with consistent accuracy across the evaluated physical quantities. The agreement in liquid fraction and temperature confirms the ability of the model to resolve both conduction-controlled and convection-enhanced melting stages. The quantitative comparison further supports the suitability of the numerical implementation for the subsequent reliability-oriented analysis.

A grid independence study was performed to evaluate the sensitivity of the numerical solution to spatial discretization. Three grid resolutions were considered: a coarse grid with m, a medium grid with m, and a fine grid with m. The comparison was based on final liquid fraction, stored energy, and energy deficit.

The grid deviation was calculated as:

where ϕ represents the evaluated numerical quantity.

Figure 6 shows the grid independence analysis, including thermal response, liquid fraction evolution, stored energy, computational cost, and the summary of grid-dependent variations.

As shown in Figure 6, the solution exhibits clear convergence with mesh refinement. The variation between the medium and fine grids remains small for the main thermal and energy-based outputs, while the computational cost increases substantially for the finest grid. The medium grid therefore provides an appropriate balance between numerical accuracy and computational efficiency.

The sensitivity of the solution to temporal discretization was assessed using three time-step sizes: s, s, and s. The comparison focused on the same key outputs used in the grid independence analysis, including final liquid fraction, stored energy, and energy deficit.

The time-step deviation was calculated as:

where denotes the result obtained using the smallest time step.

Table 5 summarizes the combined grid independence and time-step sensitivity results. The table shows that the selected medium grid and time step provide stable results with deviations below 2% relative to the finest spatial and temporal configurations.

Based on these results, the grid size m and the time step Δt = 0.1 s was selected for all production simulations.

Figure 7 presents the time-step independence analysis and illustrates the convergence behavior of the numerical solution with decreasing time-step size.

As shown in Figure 7, reducing the time step below 0.1 s produces only marginal changes in the evaluated quantities, while increasing computational demand. This confirms that the selected time step is sufficient to resolve the transient phase-change behavior and the associated energy-based outputs.

Energy conservation was verified by comparing the cumulative input energy with the sum of stored and extracted energy over the full daily cycle. The global balance error was calculated as:

The prescribed acceptance criterion was:

Table 6 presents the energy conservation results over the daily cycle. The balance error remains below 0.013 at all reported times and decreases to 0.006 at the end of the cycle.

The results confirm that the numerical model preserves the global energy budget with satisfactory accuracy. The small balance errors indicate that numerical diffusion and discretization effects do not significantly affect energy tracking.

Figure 8 provides a time-resolved representation of the energy conservation behavior, including cumulative energy balance, instantaneous energy rates, residual evolution, and summary energy quantities.

As shown in Figure 8, the cumulative input energy is consistently balanced by stored and extracted energy. The residual remains small throughout the cycle, confirming that the numerical scheme does not introduce artificial energy gain or loss. This verification is essential because the proposed reliability metrics are directly derived from energy quantities.

The numerical solution was also evaluated for physical boundedness and stability. Throughout all simulations, the liquid fraction remained within the physically admissible range . The output power remained non-negative, and the cumulative available energy increased monotonically with time. The reliability metrics also remained bounded, with and . No non-physical oscillations or numerical instabilities were observed during the simulations.

The validation and verification analyses demonstrate that the numerical model is suitable for evaluating PCM-based thermal energy storage under intermittent forcing. The benchmark comparison confirms that the model reproduces the main features of n-octadecane melting with relative errors below 4%. The grid and time-step sensitivity analyses show that the selected discretization parameters provide stable and converged results. The energy conservation check confirms that the global energy budget is preserved over the full daily cycle. Together, these results provide a reliable numerical basis for the energy-based reliability assessment presented in the following sections.

The melting process begins in a conduction-dominated regime, in which heat transfer is governed primarily by thermal diffusion in the vicinity of the heated wall. During this initial stage, the liquid fraction increases gradually and remains spatially localized. As the temperature field evolves, buoyancy effects become increasingly important, leading to the development of natural convection currents within the molten region. This transition marks a shift to convection-dominated heat transfer and is associated with a substantial increase in the melting rate.

Although the same qualitative sequence is observed across all forcing scenarios, the timing and intensity of the transition depend strongly on the level of input intermittency. Increasing fluctuation intensity delays the onset of convection and reduces the effectiveness of heat transport within the enclosure.

At the end of the cycle, the final liquid fraction decreases from 0.933 in Scenario I to 0.921 in Scenario II and 0.904 in Scenario III, indicating a progressive reduction in phase-change extent. The scenario-dependent timing of the transition to convection-dominated melting is quantified in Section 5.6.

Figure 9 shows the evolution of the phase-change process under different forcing scenarios.

As shown in Figure 9, the melting process follows a similar qualitative trajectory in all cases. However, increasing intermittency leads to a delayed development of convective flow and a reduced overall melting rate. This behavior indicates that temporal variability limits the effective utilization of the available thermal energy by constraining internal heat transport.

The total daily input energy remains fixed at 103.13 kJ across all scenarios, while the effective storage capacity lies in the range of 52–57 kJ, confirming that the system operates in an input-surplus regime. Despite this apparent sufficiency, the ratio of available energy to daily demand decreases systematically from 0.958 in Scenario I to 0.914 in Scenario II and 0.856 in Scenario III.

This result demonstrates that average input energy alone does not guarantee reliable performance. The temporal distribution of energy input plays a critical role in determining how much of the stored energy can be effectively delivered to the load.

Figure 10 presents a time-resolved comparison between available energy and demand.

As shown in Figure 10, the input power and demand profiles are not temporally aligned, leading to alternating periods of energy surplus and deficit. During surplus periods, energy is accumulated within the PCM, while during deficit periods, the system relies on previously stored energy. However, the cumulative energy curves reveal that temporal mismatch prevents full demand satisfaction, even when the total input energy is sufficient.

The principal performance metrics for all scenarios are summarized in Table 7.

The results reveal a consistent degradation in all reliability-oriented performance metrics as the level of renewable intermittency increases. While the final liquid fraction decreases only slightly from 0.933 to 0.904, substantially larger changes are observed in demand-oriented indicators. ERI decreases by approximately 10.5% and LMI decreases by approximately 14.0%, whereas the energy deficit increases more than sevenfold and the critical duration of unmet demand nearly triples. These findings indicate that thermal indicators alone may underestimate the operational impact of intermittency and highlight the importance of evaluating energy storage systems using reliability-based performance measures.

The reliability metrics exhibit a consistent degradation as intermittency increases. The Energy Reliability Index decreases from 0.962 in Scenario I to 0.917 in Scenario II and 0.861 in Scenario III. A similar trend is observed for the Load Matching Index, which decreases from 0.931 to 0.874 and 0.801.

Figure 11 illustrates the variation of the reliability metrics across the three forcing scenarios.

The figure compares the response of the Energy Reliability Index (ERI), Load Matching Index (LMI), cumulative energy deficit, and duration of unmet demand across the three forcing scenarios.

Table 8 summarizes the ensemble statistics of the four reliability-oriented performance metrics used throughout the study. The reported means, standard deviations, and 95% confidence intervals demonstrate the statistical stability of the results and confirm that the observed differences among the intermittency scenarios are robust. The statistical significance of the differences between scenarios is evaluated in Table 9.

The consistently low p-values and large effect sizes confirm that the observed degradation in reliability metrics is statistically significant and reflects systematic changes in system behavior rather than sampling variability.

The energy deficit increases from 0.82 kJ in Scenario I to 2.74 kJ in Scenario II and 5.92 kJ in Scenario III. Similarly, the critical duration of unmet demand increases from 2.1 h to 3.8 h and 5.6 h.

Figure 12 presents the temporal characteristics of energy deficit under different forcing conditions.

As shown in Figure 12, the deficit is characterized by both magnitude and duration. Increasing intermittency results in longer consecutive deficit periods and larger cumulative deficits. This behavior is particularly pronounced during peak demand intervals, where insufficient energy input leads to sustained unmet demand.

Intermittency introduces a structural modification in system behavior. Compared with the deterministic case, the onset of convection is delayed by 0.8 h in Scenario II and 2.1 h in Scenario III. At the same time, the peak output power decreases from 242 W to 226 W and 204 W, respectively.

These changes are accompanied by a reduction in the available-energy-to-demand ratio from 0.958 to 0.914 and 0.856. Together, these results indicate that temporal fluctuations affect both the rate of thermal development and the system’s ability to meet demand.

The observed degradation in performance can be explained by the interaction between phase-change dynamics and the temporal structure of energy input. During low-input periods, the imposed heat flux is insufficient to sustain buoyancy-driven flow, and heat transfer remains conduction-dominated. This limits the growth of the liquid region and reduces the accumulation of usable thermal energy.

During high-input periods, the system is unable to fully utilize the incoming energy due to thermal inertia and limited convective transport. As a result, part of the input energy becomes effectively inaccessible from a system-level perspective.

The onset of convection is identified using a velocity-based criterion, defined as the time at which the spatially averaged velocity magnitude exceeds m/s for at least 300 s. The corresponding transition times are reported in Table 10.

The results show a systematic delay in the transition to convection with increasing intermittency. This delay reduces the duration of the convection-dominated regime, leading to lower heat transfer rates and reduced energy availability.

To isolate the contribution of the proposed framework, three evaluation layers are compared: thermal-only analysis, energy-based analysis, and the full reliability-oriented framework.

The results show that thermal indicators alone produce only minor differences between scenarios, whereas energy-based and reliability-based metrics reveal substantially larger performance separation. This confirms that conventional thermal analysis is insufficient for capturing the operational impact of intermittency.

The robustness of the proposed framework is evaluated using alternative demand profiles. The results are summarized in Table 12.

The results indicate that the qualitative trends remain consistent across different demand schedules, confirming that the reliability metrics are not sensitive to a specific demand configuration.

The results demonstrate that optimizing PCM-based thermal energy storage systems requires consideration of temporal characteristics in addition to traditional thermal performance criteria. System design should account for the interaction between input variability, demand profiles, and internal heat-transfer dynamics.

Figure 13 provides a synthesized interpretation of system behavior across different intermittency levels.

As shown in Figure 13, the system transitions from a reliability-dominated regime at low intermittency to a deficit-dominated regime at high intermittency. This transition highlights the importance of considering temporal variability in the design and evaluation of thermal energy storage systems. These findings also indicate that reliability-oriented PCM storage models could be coupled with advanced building-energy control frameworks to support climate-resilient operation under thermal stress and variable renewable supply (Samaei and Riffat 2026b).

The comparison in Table 13 shows that recent PCM-TES studies have mainly evaluated performance using thermal, energy, exergy, or material-based indicators. These indicators are valuable for characterizing heat-transfer behavior and storage capacity, but they do not directly quantify whether the stored energy is available when demand occurs. The present study differs by introducing demand-oriented reliability metrics into a physics-resolved PCM model. This enables the effect of renewable intermittency to be quantified not only through changes in liquid fraction or stored energy, but also through demand satisfaction, temporal load matching, cumulative energy deficit, and critical duration of unmet demand. This distinction is important because the results show that final liquid fraction changes only slightly across the intermittency scenarios, whereas reliability indicators reveal a much stronger degradation in operational performance.

This study developed a reliability-oriented framework for evaluating PCM-based thermal energy storage systems under intermittent renewable forcing. The proposed approach combines a physics-resolved phase-change model with an energy-based performance layer, enabling direct assessment of the ability of thermal energy storage systems to satisfy time-dependent demand over a complete operating cycle.

A key contribution of the study is the transition from conventional temperature- and heat-transfer-based evaluation toward demand-oriented energy assessment. The concept of available energy provides a physically consistent measure of the portion of stored thermal energy that can be effectively delivered to the load. Based on this concept, the Energy Reliability Index (ERI) and Load Matching Index (LMI) quantify demand satisfaction and temporal supply-demand alignment, while cumulative energy deficit and critical duration characterize the magnitude and persistence of unmet demand.

The results demonstrate that systems subjected to identical total input energy of 103.13 kJ can exhibit substantially different levels of performance under varying degrees of renewable intermittency. Increasing the fluctuation intensity from σ = 0.00 to σ = 0.40 reduces ERI from 0.962 to 0.861 and LMI from 0.931 to 0.801. Simultaneously, the cumulative energy deficit increases from 0.82 kJ to 5.92 kJ, while the critical duration of unmet demand extends from 2.1 h to 5.6 h. These changes occur despite relatively small variations in final liquid fraction, indicating that conventional thermal indicators alone are insufficient for evaluating operational performance under intermittent renewable input.

The analysis further reveals that intermittency alters the internal phase-change dynamics by delaying the transition from conduction-dominated to convection-dominated melting and reducing the effectiveness of convective heat transfer. Consequently, a portion of the incoming energy becomes unavailable during critical demand periods because of temporal mismatch between energy input, storage processes, and demand requirements.

The findings demonstrate that system performance is governed not only by the total amount of stored energy but also by the timing at which that energy becomes available relative to demand. Demand satisfaction decreases as intermittency increases, while both the magnitude and duration of unmet demand grow significantly. These results highlight the importance of incorporating temporal variability into the design, assessment, and optimization of PCM-based thermal energy storage systems intended for renewable-integrated energy applications.

The present study considers a two-dimensional PCM enclosure with a single PCM material and prescribed demand profiles. Although this configuration enables systematic investigation of intermittency effects, practical thermal energy storage systems may involve more complex geometries, multiple PCM layers, cascaded storage arrangements, and variable operating conditions. Future research should extend the proposed framework to three-dimensional systems, different PCM materials, hybrid storage configurations, and experimentally validated renewable-demand datasets. In addition, integration with techno-economic analysis and long-term operational assessment would provide further insight into the practical deployment of reliability-oriented thermal energy storage systems.