A novel two-dimensional effectiveness-NTU method for high-temperature cascade latent heat storage

Research Paper A novel two-dimensional effectiveness-NTU method for high-temperature cascade latent heat storage Hussein Alawai Ibrahim Al-Saaidi , Anton Lopez-Roman , Cristina Prieto * University of Seville, Energy Engineering Department, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain ARTICLE INFO Keywords: Thermal energy storage system Cascade phase change materials Twodimensional ε -NTU method CFD model Metal wool ABSTRACT This study presents a novel two-dimensional effectiveness number of transfer unit ( ε -NTU) method to characterise and optimise cascade latent heat storage (CLHS) systems, incorporating metal wool enhancement for improved thermal performance. The proposed model provides a more accurate and computationally efficient alternative to conventional one-dimensional approaches for analysing heat transfer in high-temperature latent heat storage (LHS) applications. Current performance evaluation methods still lack precision due to the reliance on simplified assumptions regarding material properties and heat transfer dynamics. Moreover, most enhancement strategies focus on single-PCM systems, necessitating dedicated research efforts for CLHS-specific solutions. Validation through computational fluid dynamics (CFD) simulations demonstrates strong agreement, confirming the method’s reliability in predicting thermal behaviour during both charging and discharging processes. The results reveal that the cascade configuration reduces thermal resistance and enhances heat transfer efficiency, while the addition of metal wool increases the effectiveness of the system by up to 54% in single latent heat storage (SLHS) and 20% in CLHS without additives. These findings establish the two-dimensional ε -NTU method as a robust design and optimisation tool for next-generation thermal energy storage systems in renewable energy real-world applications. 1. Introduction Solar energy is one of the most promising renewable energy resources. However, its intermittent nature poses a significant challenge to ensure a stable energy supply. To address this problem, thermal energy storage (TES) systems have been developed as an effective solution to improve energy efficiency and balance supply and demand fluctuations. TES systems are increasingly being used to store excess energy generated from renewable sources, such as solar or wind, and release it when demand is high, thus improving energy efficiency and grid stability [1,2]. TES can be classified into three main categories: sensible heat storage, latent heat storage (LHS), and thermochemical heat storage. Among these, LHS systems stand out due to their high energy storage density and ability to maintain a constant phase change temperature [3,4]. Despite these advantages, a major limitation of LHS systems is the low thermal conductivity of phase change materials (PCMs), which hinders heat transfer and leads to a slow charging and discharging rate [5,6]. To mitigate these drawbacks, cascaded latent heat storage (CLHS) has emerged as a promising alternative, offering improved heat transfer rates, enhanced efficiency, and a more uniform outlet temperature of the heat transfer fluid (HTF) compared to single PCM systems [7]. However, CLHS technology remains largely in the theoretical and laboratory research stages, with limited real-world applications [8]. Several studies have explored different CLHS configurations and their thermal performance enhancements. Farid and Kanzawa [9] investigated a thermal storage system with three PCMs, reporting a 15 % increase in heat transfer rate compared to single PCM systems. Similarly, Michels and Pitz-Paal [10] conducted experimental and numerical assessments of CLHS for high-temperature applications, demonstrating its benefits in achieving higher heat transfer rates and more stable outlet temperatures. However, they highlighted that the design complexity and low thermal conductivity of PCMs remain significant obstacles. Further numerical studies have explored innovative CLHS configurations. Rudra Murthy et al. [11] compared the heat transfer performance of a tapered CLHS shell-and-tube system with a conventional cylindrical design, reporting a 17.6 % improvement in mean power output and a higher melting rate. Likewise, Elsanusi and Nsofor [12] studied the thermal behaviour of multiple PCMs in a heat exchanger and found that a series PCM arrangement reduced the total melting time by 15.5 % compared to a single PCM setup. For the design and optimisation of thermal energy storage systems, * Corresponding author. E-mail address: [email protected] (C. Prieto). Contents lists available at ScienceDirect Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng https://doi.org/10.1016/j.applthermaleng.2025.127452 Received 22 April 2025; Received in revised form 30 June 2025; Accepted 3 July 2025 Applied Thermal Engineering 278 (2025) 127452 Available online 5 July 2025 1359-4311/© 2025 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). the one-dimensional efficiency-number of transfer units ( ε -NTU) method is widely used. Tay et al. [13] introduced a one-dimensional ε -NTU technique to characterise PCM-based TES systems, allowing rapid prediction of the heat transfer rate and optimisation of the system. Furthermore, [14] validated this method using a CFD-based shell-andtube PCM model, predicting both the melting and solidification processes. Despite advancements in the one-dimensional ε -NTU method, this method assumes uniform thermophysical properties and simplified heat transfer behaviour [15,16] and which makes the selection of PCMs in CLHS systems more complicated due to uneven melting and solidification time occurring due to the different thermophysical properties for these materials. These approximations introduce significant errors in real-world applications and limit their applicability to more complex CLHS configurations, highlighting the need for more advanced modelling techniques. There are few theoretical studies that thoroughly examine the optimisation and performance of the CLHS system. Xu and Zhao [17] optimised the CLHS design using entransy theory, demonstrating that CLHS configurations could extend the applicable temperature range for multigrade thermal energy applications. Furthermore, in [18] both thermodynamic irreversibility and heat transfer rate methods based on entropy and entransy theories were used to optimise the CLHS system. The results indicated that thermal efficiency in entransy optimisation surpasses that in entropy optimisation, whereas exergy efficiency is higher in entropy optimisation compared to entransy optimisation. Nevertheless, a challenge with entransy dissipation theory lies in the incomplete definition of stored and transferred heat, as noted by Kostic et al. [19]. Metal wool has been proposed as an effective heat transfer enhancement strategy as a result of its high thermal conductivity, ease of integration with PCMs, and cost-effectiveness. Prieto et al. [20] demonstrated that incorporation of metal wool into PCM systems increased the effective thermal conductivity by 300 %, significantly improving heat transfer performance in high-temperature applications. Similarly, Favache et al. [21] investigated the use of high-conductivity metallic wool as a heat transfer augmentation technique in PCMs, reporting improved heat transfer during solid-phase transitions, although natural convection effects in the liquid phase remained limited. Prieto et al. [22] explore the feasibility of employing metal wool as an economical method to improve thermal conductivity for latent thermal energy storage (TES) in solar process heat applications. To ensure the stability of metal wool in typical industrial conditions, experiments were conducted under high temperatures and in an inert atmosphere. These tests aimed to assess the thermal degradation of metal wool when subjected to elevated thermal cycling temperatures. The tests were carried out at a broad range of temperatures 200 C to 500 C in an inert atmosphere. The results reveal no chemical or physical degradation, thus confirming the suitability of this technique for the considered application. Despite these promising findings, key research gaps remain in the optimisation of the CLHS system. Current performance evaluation methods still lack precision due to the reliance on simplified assumptions regarding material properties and heat transfer dynamics. Moreover, most enhancement strategies focus on single-PCM systems, necessitating dedicated research efforts for CLHS-specific solutions. In this paper, we propose a two-dimensional ε -NTU technique to address these challenges, offering a more accurate and computationally efficient approach to analyse high-temperature CLHS systems. Unlike conventional one-dimensional methods, this technique accounts for spatial variations in heat transfer and material properties, which makes it better suited for complex TES applications. To validate this approach, a computational fluid dynamics (CFD) model is implemented, allowing for precise comparison and performance assessment. Furthermore, we evaluated the thermal performance of metal wool-enhanced CLHS configurations, examining their impact on heat transfer efficiency, charging/discharging rates, and overall system effectiveness. The findings of this study contribute to the development of more robust and Nomenclature A area of heat transfer, m 2 Cp PCM specific heat of the PCM C p HTF specific heat HTF, kJ/kg K SLHS single latent heat storage CLHS cascaded latent heat storage h f heat transfer coefficient of the HTF, W/m 2 K h sensible enthalpy, kJ/kg H total enthalpy, kJ/kg k PCM thermal conductivity of the PCM, W/m K k PCM1 thermal conductivity of the PCM 1 , W/m K k PCM2 thermal conductivity of the PCM 2 , W/m K k w thermal conductivity of the tube wall, W/m K k e effective thermal conductivity, W/m K L total length of the tube, m L 1 length of the first container includes PCM 1 , m L 2 length of the second container includes PCM 2 , m LHS latent heat energy storage system − NTU number of transfer unit N number of PCMs P P factor, − Q act actual stored energy within a PCM system, kJ Q max maximum energy storage, kJ r i inner radius of the tube, m r o outer radius of the tube, m r max radius of PCMs, m R 1 total thermal resistance of the first parallel heat flow path of HTF, tube, and PCM 1 , K/W R 2 total thermal resistance of the second parallel heat flow path of HTF, tube and PCM 2 , K/W R PCM1/PCM2/PCMn thermal resistance of the isothermal heat flow path consisting of PCM 1 , PCM 2 and PCM n , K/W R HTF HTF thermal resistance, K/W R HTF1 HTF thermal resistance for the first heat flow path, K/W R HTF2 HTF thermal resistance for the second heat flow path, K/W R iso total thermal resistance using isothermal heat flow, K/W R parallel total thermal resistance using parallel heat flow, K/W R PCM PCM thermal resistance, K/W R PCM1 PCM thermal resistance for the first heat flow path, K/W R PCM2 PCM thermal resistance for the second heat flow path, K/W R T total thermal resistance, K/W R Tube tube thermal resistance, K/W R Tube1 tube thermal resistance for the first heat flow path, K/W R Tube2 tube thermal resistance for the second heat flow path, K/W R Tuben tube thermal resistance for the number of heat flow path, K/W r Stefan–Boltzmann constant, W/m 2 K 4 γphase change fraction, − ε heat exchanger effectiveness, − ε s emissivity of the radiating surface βcoefficient of thermal expansion, 1/K U overall heat transfer coefficient, W/m 2 K ˙ mmass flow rate of HTF, kg/s δporosity, − λthe ratio of the intersection radius to diameter of the wool fibre H.A. Ibrahim Al-Saaidi et al. Applied Thermal Engineering 278 (2025) 127452 2 scalable CLHS designs, which ultimately supports the advancement of next-generation thermal energy storage solutions for renewable energy applications. 2. Mathematical formulation A mathematical model based on one dimensional ε -NTU technique has been developed and experimentally validated for the tube in a PCM system using a computational fluid dynamics (CFD) model by Tay et al. [13]. The one-dimensional method cannot be implemented to design complicated configurations such as a CLHS and LHS with enhancement techniques because it assumes to simplify the mathematical model by making the heat flow only in one direction. Moreover, the twodimensional method can be a more appropriate representation of the thermal resistance of a CLHS that includes different thermophysical properties and volumes of PCMs. This approach helps to determine the actual average heat transfer rate and the phase-transition time throughout the phase-transition process, using a specific set of design parameters for the tube-in-tank configuration. Average effectiveness was described as the heat flow to a heat source/sink of infinite specific heat and was expressed as the following: ε =1−exp(− NTU) = Qact Qmax (1) The NTU at any point in time can be defined as [23]: Fig. 1. (a) Parallel and (b) isothermal heat flow of the CLHS system. Fig. 2. (a) Parallel heat flow model, and (b) thermal circuit. H.A. Ibrahim Al-Saaidi et al. Applied Thermal Engineering 278 (2025) 127452 3 NTU =UA (˙ mCp)=1 (RT˙ mCp)(2) The ε -NTU method depends on the characterisation of the thermal resistance between the HTF and the phase change boundary. An accurate depiction of the thermal resistance in a CLHS system would aid in the design and optimisation of such a system without requiring CFD modelling. Consequently, it is suggested to utilise CFD modelling to calculate the average effectiveness, which will then be used to calculate the equivalent total thermal resistance of the CLHS system and from this value the equivalent total thermal resistance of the CLHS system. Two-dimensional heat transfer is bounded by the limits of parallel heat flow, characterised by an infinite transverse thermal resistance, and isothermal heat flow, where the transverse resistance is zero. For a CLHS, parallel heat flow occurs when there is no transverse heat flow, thus allowing heat flow only in a one-dimensional direction parallel to the tube wall, as shown in Fig. 1(a). However, isothermal heat transfer, characterised by the absence of lateral thermal resistance, allows heat flow in a one-dimensional direction perpendicular to the tube wall, as shown in Fig. 1(b). The actual phase change process is a combination of these two mechanisms. Gorgoleski [24] developed the mathematical model for the concept of steel frames-bridged insulation. To determine total thermal resistance, a factor called P is used to represent the proportion of resistance on parallel and isothermal paths, with values ranging from 0 to 1. Typically, the factor P is assumed to be 0.5 according to Belusko [25], but it should ideally be determined through experimentation or threedimensional conduction modelling, such as computational fluid dynamics (CFD). Actual thermal resistance, R T , is identified by Eq. (3) [26]. It is suggested to employ CFD to assess the total thermal resistance and then identify the suitable P factor. RT=P⋅Rparallel +(1−P)⋅Riso (3) The total thermal resistance of the system, when the heat flow paths are run parallel, is referred to as the parallel resistance, as explained in Eq. (4). Fig. 2 illustrates the resistance diagram that clarifies the identification of two heat flow paths. The first path involves the HTF, tube wall, and PCM 1 , while the second path consists of the HTF, tube wall, and PCM 2 . Rparallel =R1⋅R2.Rn/(R1+R2+Rn)(4) R 1 and R 2 are the thermal resistances through the first and second heat flow paths, while R n represents the thermal resistance for an unknown number of PCMs in the cascade system, as shown in the equations. (5–7), respectively. R1=RHTF1 +RTube1 +RPCM1 (5) R2=RHTF2 +RTube2 +RPCM2 (6) Rn=RHTF n +RTube n +RPCM n (7) R HTF1 and R HTF2 in Eqs. (8)–(9), denote the HTF thermal resistances for the first and second heat flow paths, while R HTFn in Eq. (10) represent the thermal resistances of the unknown number of HTFs, which is characterised by forced internal convection. RHTF1 =1/(2 π riL1hf)(8) RHTF2 =1/(2 π riL2hf)(9) RHTFn=1/(2 π riLnhf)(10) R Tube1 and R Tube2 in Eqs. (11–12) denote the thermal resistances of the tube wall for the first and second heat flow paths, while and R Tuben in Eq. (13) represents the resistance of the unknown number of tubes. RTube1 =ln(ro/ri)/(2 π kwL1)(11) RTube2 =ln(ro/ri)/(2 π kwL2)(12) RTuben=ln(ro/ri)/(2 π kwLn)(13) R PCM1 and R PCM2 in Eqs. (14)–(15) denote the PCMs resistances in the first and second heat flow paths, while R PCMn in Eq. (16) represents the unknown number of PCMs in the cascade system, defined by conduction and the relevant shape factor. RPCM1 =ln({γ(r2 max −r2 o)+r2 o}1/2/ro)/2 π L1kPCM1 (14) RPCM2 =ln({γ(r2 max −r2 o)+r2 o}1/2/ro)/2 π L2kPCM2 (15) RPCMn=ln({γ(r2 max −r2 o)+r2 o}1/2/ro)/2 π LnkPCMn(16) Fig. 3. Schematic of the physical model of (a) single (b) cascade latent heat storage. H.A. Ibrahim Al-Saaidi et al. Applied Thermal Engineering 278 (2025) 127452 4 In Equation (17), the second resistance, Riso, denotes the total thermal resistance of the system, achieved by assuming that each component of the system (HTF, tube wall, and PCM) maintains an isothermal state. This resistance is derived from the resistance circuit illustrated in Fig. 2, which indicates the presence of three distinct layers. Riso =RHTF +RTube +RPCM1 /PCM2/PCMn (17) The resistance of the first layer (HTF) is given in Eq. (18). RHTF1 =1/(2 π riLhf)(18) The resistance of the second layer (tube wall) is given in Eq. (19). RTube =ln(ro/ri)/2 π kwL(19) The resistance of the third layer (PCM 1 , PCM 2 , and PCM n ) is given in Eq. (20). This layer consists of the total number of PCMs. The thermal resistance of this section is defined by three resistances in parallel (Fig. 2 b), as revealed in the equations. (20) and (21). RPCM1 /PCM2/PCMn =1/⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ [2 π L1kPCM1/ln({γ(r2 max −r2 o)+r2 o}1/2/ro)]+ [2 π L2kPCM2/ln({γ(r2 max −r2 o)+r2 o}1/2/ro)] [2 π LnkPCMn/ln({γ(r2 max −r2 o)+r2 o}1/2/ro)] +⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (20) 1/Rpcm1 /PCM2/PCMn =1/RPCM1 +1/RPCM2 +1/RPCMn (21) 3. Problem statement and numerical modelling procedure 3.1. Physical model A three-dimensional schematic diagram for two vertical cylindrical shell and tube arrangements representing SLHS and CLHS with 90 % porosity of tungsten stainless steel metal wool is shown in Fig. 3 (a) and (b). The metal pipe is used with an inner diameter of 22.22 mm and a wall thickness of 1.65 mm, while the shell, which encloses the phase change materials (PCMs), is constructed from stainless steel with an outer diameter of 73 mm and a wall thickness of 2.11 mm. The cylindrical latent heat storage unit has a length of 888 mm. In the case of the CLHS configuration, the storage shell is divided into two sections of equal lengths, as mentioned by Jain et al. [27]. Therminol 66 is used as HTF, which passes through the inner tube from top to bottom with a uniform inlet velocity range from 0.01 to 0.12 m/s, while the inlet temperature of 615 K and 519 K is implemented during the melting and solidification processes for all simulations. The initial temperature for all simulations is set at 519 and 615 K for the melting and solidification processes, respectively. HTF discharges from the storage unit at atmospheric pressure. The adiabatic boundary conditions are implemented on the faces between stages and on the external faces that are exposed to the environment. Table 1 lists the thermophysical characteristics of PCM 1 , PCM 2 , and HTF employed in this study. 3.2. Governing equations and assumptions The enthalpy-porosity method was implemented to perform a threedimensional transient simulation of the melting and solidification processes, in which the interfacial phase change zone (mushy zone) is treated as a porous medium [28]. The boundary condition used in this study is summarised in Table 2. The equations covered in the numerical analysis are noted by Mayeli et al. 2021 [29]: Heat transfer fluid domain The continuity equation for the HTF domain is as follows: ∂ρ f ∂ t+∇•( ρ fV →)=0 (22) The momentum conservation equation for the HTF domain as follows: ρ f(DV → Dt )= − ∇p−2 3∇( μ f∇• V →)+∇•[ μ f(∇V →+(∇V →)T)] (23) The energy conservation equation for the HTF domain is as follows: ρ fcpf(DT Dt )= ∇•(kf∇T)(24) Phase change materials domain Continuity equation ∇.V →=0 (25) Where v is the velocity vector, and its components, u, v, and w, are located, respectively, in the r, Θ and z directions Momentum equation ∂ V → ∂ t+V →⋅∇V →=1 ρ (− ∇P+ μ ∇2V →+ ρ g →β(T−Tref))+Sm(26) Energy equation ∂ h ∂ t+ ∂ H ∂ t+∇⋅(V →h)= ∇⋅(k ρ Cp∇h)(27) Where, h is sensible enthalpy and H is total enthalpy and P, ρ , V and g denote the fluid pressure, density, velocity, and acceleration due to gravity, respectively. The variation of density is presented as the Boussinesq assumption: ρ = ρ l /(β(T −T l ) +1) where ρ l represents the density of liquids and β is the thermal expansion. Furthermore, S m is a source of momentum that adds to the momentum equation and is defined by Olfian et al. 2020 [30] Table 1 Thermophysical properties of PCMs, HTF [27], and metal wool [20]. Property PCMs HTF Metal wool PCM 1 (NaNO 3 ) PCM2 (NaNO 2 ) Therminol 66 Tungsten steel Temperature of melting/ solidification (K) 579 555 − − Density [kg/m 3 ] 1908 1812 775 7913 Latent heat of fusion [kJ/kg] 176 180.12 161 − Specific heat solid [kJ/ kg K] 1.60 1.733 −0.488 Specific heat liquid [kJ/kg K] 1.655 2.553 2.730  Thermal conductivity solid [W/m K] 0.8 0.765 −66 Thermal conductivity liquid [W/m K] 0.68 0.665 0.0893 − Dynamic viscosity [mPa s] 2.69 2.666 0.336 − Coefficient of Thermal Expansion (1/K) 0.0004 0.00028 − Table 2 Boundary conditions of all cases for both charge and discharge processes. Boundary conditions Melting Solidification Velocity of HTF, [m/s] 0.01–0.12 0.01–0.12 HTF inlet temperature, [K] 615 519 PCM initial temperature, [K] 519 615 H.A. Ibrahim Al-Saaidi et al. Applied Thermal Engineering 278 (2025) 127452 5 Sm= − Am(1−γ)2 (γ3+Φ)(v−vp)(28) Equation (28) establishes a distinction between the liquid and solid phases of PCM based on the enthalpy-porosity technique, Gowreesunker et al. [31]. A m , the mushy zone constant, with a value between 10 4 and 10 7 , υ p is the solid velocity, and Φ is a small number of 0.001 that prevents zero division. γ is the liquid fraction that results from the transition between the liquid and solid phase, which is used to characterisee the degree of melting process is a constant between 0 and 1 and can be described as follows: γ=⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 0 T <Tsolidus T−Tsolidus Tliquidus −Ts Tsolidus ≤T≤Tliquidus 1 T >Tliquidus (29) The results from the solver are then analysed in the post-processing stage using the same CFD programme to give the results. The theoretical model proposed for the effective thermal conductivity of composites, introduced by [32] and experimentally validated in [33], is used as follows: ke= 2 √ Ra+Rb (30) Ra=4λξ kf+[2 π λ2ξ2 3+ π ξ2 2(1 2λξ −3+2 2 √)](ks−kf) (31) Rb= 2 √−4λξ kf+ 2 √ π ξ2(ks−kf)(32) ξ=54CA2−4B3+6 3 √[(27C2A2−4CB3)1 2A]1 3 3A + 4B2{54CA2−4B3+6 3 √[(27C2A2−4CB3)1 2A](−1 3)}−B 3A (33) A=4 2 √ π λ3 3−3 2 √ π λ,B=3 π  2 √,C=1−δ(34) The effective thermal conductivity of the composite was expressed through the conductivities of the fluid (kf) representing PCM and the thermal conductivity of the solid structure (ks) representing metal wool, the porosity of the metal wool (δ) and the ratio of the radius intersection to the ligament diameter of wool fiber (λ =r/d).).). This parameter was derived from experimental observations. According to the data in [33], the value λ =60 provided the best fit. The following presumption was considered when analysing processes using the LHS system: •Except for the density, all the physical properties of PCM are constant. •The HTF characteristic properties are constant with temperature. •Liquid PCM is considered as a laminar, Newtonian, and incompressible fluid. •Ignoring the viscous dissipation of energy. •In both processes, it is assumed that the system is adiabatic. •The liquid phase through the PCM is modelled using the Boussinesq approximation [34]. 3.3. Numerical model The computational process can be classified into three steps: preprocessing, solver (processing), and post-processing. The pre-processing stage involves constructing the computational domain that contains the assembly of the model, meshing, and boundary conditions. To ensure accurate results with a reasonable computing cost, an independent study of the grid and the time steps needs to be conducted. In the current study, three different numbers of elements (196500, 697500, 813750) and three different time steps (0.1, 0.5, and 1 s) were evaluated. Based on analysis, number of element and time steps 697,500 and 1 are taken on the considered in the further computations as seen in Fig. 4 (a) and (b). A grid with N =697,500 cells was selected for its optimal computational cost saving. Fig. 4 (a) indicates that refining this grid further does not significantly affect the numerical solution. Additionally, as shown in Fig. 4 (b), the time step of 1 s was determined to be sufficient to maintain the solution’s independence and stability throughout the entire phase transition. The solver stage involves exporting the model to ANSYS FLUENT 2023R2 for setup and analysis. Due to the time dependence problem, a transient state was selected. An implicit secondorder transient formulation is applied to the solution. The energy had convergence criteria of 10 -9 , while momentum in all directions had a criterion of 10 -6 . According to numerous running tests, the set of underrelaxation parameters is suitable in this study to prevent divergences. The values for momentum, energy, and liquid fraction were 0.3, 0.9, and 0.1, while those for body forces, density, and pressure were 1, 0.8, and 0.3. The pressure–velocity coupling and the Semi-Implicit Method for pressure-linked equations (SIMPLE) technique are used. Pressure Fig. 4. Effect of (a) the number of computational elements and (b) time step on the evaluation of liquid fraction for NaNO 2 within computational domain. H.A. Ibrahim Al-Saaidi et al. Applied Thermal Engineering 278 (2025) 127452 6 interpolation employed the Staggering Pressure Option (PRESTO) methodology pressure correction [35]. 4. Results and discussion 4.1. Validation of CFD model To verify the accuracy of the current numerical model in predicting the phase transition process, the developed model is compared with the numerical result for SLHS and CLHS of Jain et al. [27]. They performed numerical investigations on the NaNO 3 SLHS based single stage storage with the same dimensions mentioned in detail [27], which has a shell made of steel (SS316) and HTF tube made of copper. The length of the storage is considered 888 mm. The NaNO 3 /NaNO 2 CLHS used the same geometry, but it is divided into two equal parts to establish two stages configuration. The same thermophysical properties and initial conditions were used for validation such as inlet temperature and velocity of HTF were 615 K and 0.04 s/m, respectively, as well as the initial temperature of PCMs was 519 K. Grid, time step independence and solver setting are chosen as the same reference mentioned above. Fig. 5 (a) and (b) show the comparison of progress in liquid fraction over time and temperature contours for two periods 120 mins and 240 mins between the present model and the results by [27] for SLHS and CLHS −based PCMs. The figure clearly shows that the present model is capable of accurately predicting the evaluation phase change and temperature profile. 4.2. Validation of two-dimensional ε -NTU method Figs. 6 and 7 show the comparison between the average effectiveness across different mass flow rates to area from the developed twodimensional ε -NTU method and the CFD model during the charging Fig. 5. Comparison of (a) Liquid fraction and (b) Temperature contours with results of Jain et al. [27]. Fig. 6. Validation of the melting results of ε -NTU method with CFD for single and cascade PCMs. Fig. 7. Validation of the solidification results of ε -NTU method with CFD for single and cascade PCMs. H.A. Ibrahim Al-Saaidi et al. Applied Thermal Engineering 278 (2025) 127452 7 and discharging processes for PCM configurations based on SLHS and CLHS. It can be observed that the cascade configuration dramatically increases the average effectiveness during the charging process by 37% and discharging process by 34% compared with single storage and that in agreement with previous studies [36,37]. The two-dimensional ε -NTU method and the CFD model showed good agreement across all flow rates. Furthermore, this agreement was achieved for both SLHS and CLHS configurations. One of the most important limitations of the onedimensional ε -NTU method is that it does not take into account the effect of natural convection. This is due to only considering the heat flow passing through one direction and ignoring the heat flow in other directions. Moreover, the natural convection is dependent on the temperature difference in the system, while this method is intended to be temperature independent. Meanwhile the two-dimensional method can consider the effect of natural convection by determining the heat flow in two directions from comparing the results obtained from CFD. In this way, the effect of natural convection can be adjusted through the value of P factor. Therefore, initially aiming to validate the two-dimensional methods with CFD results, in this section the buoyancy effect has been ignored in the CFD model to reduce the average error as much as possible. In addition, an important limitation of the two-dimensional method is that it can take into account only some thermophysical properties such as thermal conductivity while it ignores other properties such as enthalpy. An average error was found to be approximately 3.6 % and 5% for the melting and solidification processes, respectively. Therefore, the newly developed method, which assumed phase change, has been validated to accurately describe the cascade configuration during both charging and discharging processes. 4.3. Calculated average effectiveness and thermal resistance ratio Figs. 8 and 9 illustrate the calculated average effectiveness across different mass flow rates to area for the melting and solidification processes in both single and cascade-based PCMs, with and without metal wool configurations. The results indicate that CLHS exhibits a significant improvement in thermal performance over the single LHS configuration. The addition of metal wool as an enhancement technique further improves effectiveness in both configurations during charging and discharging processes. Specifically, the average effectiveness of CLHS increases by 43 % and 44 % compared to a single PCM configuration. When metal wool is incorporated, the effectiveness is enhanced by 54 % and 54.7 % over SLHS and by 20 % and 16 % over CLHS-based PCMs for the melting and solidification processes, respectively. The improvement in heat transfer efficiency with metal wool is attributed to two key factors: (1) the high thermal conductivity of metal wool, which facilitates rapid heat transfer, and (2) the large heat transfer surface area created by the metal wool structure, enabling better thermal distribution within the PCMs. This enhancement overcomes the inherent low thermal conductivity of PCMs, allowing for more efficient phase transition processes. Figs. 10 and 11 present the average thermal resistance within the PCMs, relative to the total thermal resistance, as a function of HTF velocity for both single and cascade configurations with and without metal wool. These figures demonstrate that thermal resistance within PCMs is significantly reduced in cascade configurations, indicating that the Fig. 8. Comparison the melting results of ε -NTU method for SLHS and CLHS with and without metal wool. Fig. 9. Comparison the solidification results of ε -NTU method for SLHS and CLHS with and without metal wool. Fig. 10. Comparison the thermal resistance ratio of the melting process for SLHS and CLHS with and without metal wool. Fig. 11. Comparison of the thermal resistance ratio of the solidification process for SLHS and CLHS with and without metal wool. H.A. Ibrahim Al-Saaidi et al. Applied Thermal Engineering 278 (2025) 127452 8 dominant resistance shifts toward HTF convection resistance (R HTF ). Additionally, a substantial reduction in resistance is observed with the introduction of metal wool in both SLHS and CLHS configurations. A notable transition is observed at a HTF velocity of 0.05 m/s, where the Reynolds number reaches values characteristic of turbulent flow. This transition results in a higher Nusselt number, thereby enhancing convective heat transfer and reducing total thermal resistance. The impact of cascade resistance varies, averaging 15 % during solidification and reaching a maximum of 80 % during melting. Notably, the thermal resistance of the tube and convective resistance within the fluid remain significant contributors to overall system performance. These findings underscore the importance of cascade configurations and metal wool integration in improving thermal management and efficiency in latent heat storage systems, with direct implications for hightemperature thermal energy storage applications. 4.4. Thermal performance assessment of SLHS and CLHS-based PCMs This section evaluates the charging and discharging processes in both single and cascade latent heat storage (LHS) systems to analyse the effect of cascaded PCM configurations. Additionally, the impact of metal wool as a cost-effective enhancement technique for high-temperature applications is examined [20,38]. The thermal performance of SLHS using NaNO 3 and NaNO 2 as storage mediums, and CLHS incorporating both NaNO 3 and NaNO 2 with and without metal wool, is assessed. Both configurations are arranged in a vertical orientation, where the heat transfer fluid (HTF) flows through the inner tube while the PCMs are contained in the surrounding outer shell. In all cases, hot HTF enters the system at 0.08 m/s from the top and exits at the bottom, following the direction of gravity, as illustrated in Fig. 3. Fig. 12 (a) presents the phase change interface and temperature distribution at different time steps during the charging phase for CLHS. Initially, the phase change occurs in PCM2, as its melting temperature (555 K) is lower than the HTF temperature. During the early phase, heat conduction governs the process, causing PCM1 to undergo a sensible heating phase before reaching its higher melting temperature (579 K). Due to the lower temperature difference between PCM1 and HTF, the energy transfer rate declines, leading to an extended charging time. As charging progresses, melting occurs more rapidly at the top than at the bottom of both stages. This is attributed to HTF temperature reduction along the downward flow and buoyancy effects, which cause the heated melt to rise, slowing the melting process at the bottom. Fig. 12. Liquid fraction and temperature distributions contours during charging for both (a) CLHS and (b) CLHS +MWconfigurations throgh different time steps. H.A. Ibrahim Al-Saaidi et al. Applied Thermal Engineering 278 (2025) 127452 9