A 1D Reduced-Order Model (ROM) for a Novel Latent Thermal Energy Storage System

Phase change material (PCM)-based thermal energy storage (TES) systems are widely used for repeated intermittent heating and cooling applications. However, such systems typically face some challenges due to the low thermal conductivity and expensive encapsulation process of PCMs. The present study overcomes these challenges by proposing a lightweight, low-cost, and low thermal resistance TES system that realizes a fluid-to-PCM additively manufactured metal-polymer composite heat exchanger (HX), based on our previously developed cross-media approach. A robust and simplified, analytical-based, 1D reduced-order model (ROM) was developed to compute the TES system performance, saving computational time compared to modeling the entire TES system using PCM-related transient CFD modeling. The TES model was reduced to a segment-level model comprising a single PCM-wire cylindrical domain based on the tube-bank geometry formed by the metal fin-wires. A detailed study on the geometric behavior of the cylindrical domain and the effect of overlapped areas, where the overlapped areas represent a deviation from 1D assumption on the TES performance, was conducted. An optimum geometric range of wire-spacings and size was identified. The 1D ROM assumes 1D radial conduction inside the PCM and analytically computes latent energy stored in the single PCM-wire cylindrical domain using thermal resistance and energy conservation principles. The latent energy is then time-integrated for the entire TES, making the 1D ROM computationally efficient. The 1D ROM neglects sensible thermal capacity and is thus applicable for the low Stefan number applications in the present study. The performance parameters of the 1D ROM were then validated with a 2D axisymmetric model, typically used in the literature, using commercially available CFD tools. For validation, a parametric study of a wide range of non-dimensionalized parameters, depending on applications ranging from pulsed-power cooling to peak-load shifting for building cooling application, is included in this paper. The 1D ROM appears to correlate well with the 2D axisymmetric model to within 10%, except at some extreme ranges of a few of the non-dimensional parameters, which lead to the condition of axial conduction inside the PCM, deviating from the 1D ROM.


Introduction
With the increasing demand for efficient energy systems to reduce the mismatch between supply and demand of energy, there is a need to develop efficient thermal energy storage (TES) systems that can store energy from time-variable sources and match timevariable demands with a constant output source [1]. The interest in energy storage devices based on latent energy storage using phase change materials (PCMs) is rapidly growing due to their isothermal nature and high energy storage density [2], resulting in compact and reduced-weight, lower cost systems [3,4] when compared to sensible thermal energy storage systems. Since there are many types of PCMs, which melt at widely varying temperature ranges, they have been used in diverse applications, such as solar heating [2], peak-load shifting for building cooling [3], pulsed-power applications [4], and space energy storage applications [5].
However, the main drawbacks of PCMs are their low thermal conductivity-in the range of 0.22 to 0.33 W/m K [6,7]-sub-cooling and phase-segregation issues [8], volumeexpansion issues, and expensive encapsulation processes that sometimes comprise up to 99% of the total PCM cost [9]. The present study overcomes these challenges by utilizing a low-cost, metal-polymer composite heat exchanger (HX), consisting of metal wires embedded in the polymer structure and integrated with a PCM to yield a compact, durable, and lightweight storage device. The proposed technology leverages the results from the authors' past work on a patented cross-media approach capable of producing entirely additively manufactured, low-cost and lightweight metal-polymer composite heat exchangers [10][11][12][13][14]. The TES units based on these HXs utilize high thermal-conductivity continuous metal wires that act as fins in a tube-bank geometry, separated by lightweight and low-cost alternate polymer channels. The metal fins on the PCM side can, thus, enhance the thermal conductivity of the PCM. The expensive encapsulation process is circumvented because the PCM is separately macro-contained inside the polymer channels. In addition, as the PCMs are embedded with metal fins, the effective thermal conductivity of the PCM is expected to increase, since the high fin density of PCM-metal structures has been reported in the literature to enhance the PCM thermal conductivity up to 42 times [7]. Thus, the overall TES design parameters, such as fin density and configuration, are critical for determining the TES performance. Consequently, computational modeling should be performed to quantify these TES design parameters governing the phase-change process of the PCMs.
Several kinds of research have been reported in the literature on the numerical study of the phase-change process in general. One of the earliest approaches reported in the literature mathematically treated the phase-change problems as moving boundary problems, also called Stefan problems [15,16], and provided an exact solution [17] by solving heat equations and interface equations based on energy conservation at the phase-change boundary. However, the interface energy equations are challenging to solve analytically due to interface issues, such as non-linear nature and physical discontinuities, due to changes in properties from one phase to another. Thus, the exact solutions are valid only for simple geometries and boundary conditions. Similarly, other analytical approaches, such as calculation of variation, successive approximation, and perturbation methods [15], are also limited to simplified domains and boundary conditions. The alternative approach uses numerical techniques generally based on temperaturebased and enthalpy-based finite-difference and finite-element methods. The temperaturebased methods track the phase-change boundary by solving two separate temperature equations on the solid and liquid computational domains through energy balance at the phase change [15,18]. However, these methods are not robust and are challenging to program due to interface issues. The enthalpy-based methods don't track the phase-change interface but rather treat the problem as a typical heat conduction problem without phase change by solving a single governing equation on the combined liquid and solid computational domain. Thus, they avoid interface issues, unlike temperature-based methods. These methods are also robust and easier to program. Similar to the enthalpy-based methods, other methods commonly used are source-term-based methods and heat-capacity methods [19], velocity-transitioned modeling [20], and first and second thermodynamic law models [21]. Thus, these numerical techniques can potentially resolve the multidimensionality and complexity of versatile analytical problems, but these numerical methods are still time-consuming and computationally expensive.
The most common application of these numerical approaches is in models of energy storage processes in shell and tube storage units. Some numerical approaches use enthalpybased [22][23][24][25] methods, while others are heat-capacity-based [26][27][28]. In both approaches they formulate the PCM problem surrounding the tubes as a 2D axisymmetric problem with an inner cylinder consisting of the flow passing inside the tube. In contrast, the outer cylinder consists of the PCM. These cylindrical-based studies have been extended to model PCM cylinders enclosing finned tubes using traditional enthalpy-based approaches [29]. Some studies also use the enthalpy-based approaches directly via commercial solvers, such as Ansys Fluent [30,31]. Similar applications on such cylindrical latent heat storage systems have also been made using analytical approaches based on the second thermodynamics law [32] and moving boundary methods [30,33]. In addition to the traditionally used approaches, some other simple numerical approaches have been developed, such as the temperature and thermal resistance iteration methods [34], which are based on energy conservation and are implemented for concentric tube storage units. However, most of these transient studies are based on the 2D computational domain, which is still computationally expensive.
The present study thus resolved the above shortcomings in the literature, such as computational cost due to transient simulations based on numerical methods, complex domain, challenging programs related to analytical methods, and lack of robustness, by introducing an analytically based 1D reduced-order model (ROM). This ROM can model complex HX designs, such as the present TES design, by reducing its entire domain to a segment-level domain. The 1D ROM analytically models only at the segment level of the TES, significantly reducing the computational cost and time. This makes the ROM particularly appropriate for design optimization problems when modeling needs to be repeated many times to identify the optimum solution. Furthermore, since the segment level model of the ROM is based on 1D radial conduction inside the PCM cylinders, enveloping the metal wires in a tube-bank geometry, the 1D ROM can also be used to model other cylindrical PCM-fin-based geometries. This paper presents a geometrical analysis on the effects of the domain simplification, such as overlaps between neighboring PCM cylinders on the performance of the 1D ROM. Moreover, the 1D ROM has been found to be experimentally validated to about 15% for a specific case as in another work of the authors [35]. However, a validation of the 1D ROM for a wide range of geometric and operating parameters is still needed. The present paper discusses validation of the 1D ROM with a 2D computational fluid dynamics (CFD) axisymmetric model, commonly used in the literature for cylindrical PCM geometries for a wide range of input parameters.

Design of a Novel TES System
The novel TES system includes the cross-media-based metal-polymer composite HX (see Figure 1a), in which hot fluid (or cold fluid in the reverse heat release process) flows across alternate fluid channels (n f luid ) of width W f luid and PCM-contained polymer channels (n PCM ) of width W PCM , as shown in Figure 1c. As the hot fluid flows along the fluid channels, convective heat transfer occurs between the hot fluid and the metal wires of diameter D, exposed on the fluid side. The wires are arranged in a tube-bank geometry governed by transverse (S T ) and longitudinal (S L ) spacing parameters (see Figure 1b). The heat then continuously conducts through the wires axially, up to t wall length, until it reaches the PCM side. In this so-called cross-media approach, small-diameter continuous metal wires of high thermal conductivity, such as copper, are laid across all the fluid and PCM channels, providing a direct heat transfer route between the hot fluid and the PCM (see Figure 1c,e) [10][11][12][13][14]. Since the heat transfer route is independent of the thermal conductivity of the polymer, the resulting TES system is a low-resistance unit. When the wire reaches the PCM side, the heat starts conducting to the PCM, and when the PCM temperature reaches its phase-transition temperature (T PCM ), the PCM starts melting and stores the heat isothermally in the form of latent energy. These TES systems are entirely additively-manufactured in-house using a customized proprietary 3D printer consisting of two print heads, polymer and metal, designed and developed by a few of the co-authors. A detailed description of the fabrication technique can be found in [12,36]. two print heads, polymer and metal, designed and developed by a few of the co-authors. A detailed description of the fabrication technique can be found in [12,36]. Figure 1. Design of TES system: (a) isometric view; (b) cross-sectional view of XY plane showing metal wires in tube-bank geometry; (c) cross-sectional view of XZ plane with cross-media approach marked in dashed border rectangle; (d) picture of an additively manufactured TES system at the AHXPI laboratory [35]; (e) cross-media approach leading to melting of PCM.
The present study simulates the melting process of the PCM enveloping the metal wires in a PCM-to-fluid HX. It can be easily extended to simulate the freezing process, which is similar to the simulation of the melting process, but with minor adjustments.

1D ROM
This section includes the development of the 1D ROM for radial conduction in PCMs. It first includes the domain simplifications from the entire HX level to the segment-based level and assumptions. It then studies the condition required for the PCM to convert into a radial ROM. Lastly, it discusses the governing equations to simulate the ROM.

Domain Simplifications and Assumptions
Since the modeling of the entire HX domain is computationally expensive, it can be simplified to an independent segment-level-based model, using the following assumptions: Figure 1. Design of TES system: (a) isometric view; (b) cross-sectional view of XY plane showing metal wires in tube-bank geometry; (c) cross-sectional view of XZ plane with cross-media approach marked in dashed border rectangle; (d) picture of an additively manufactured TES system at the AHXPI laboratory [35]; (e) cross-media approach leading to melting of PCM.
The present study simulates the melting process of the PCM enveloping the metal wires in a PCM-to-fluid HX. It can be easily extended to simulate the freezing process, which is similar to the simulation of the melting process, but with minor adjustments.

1D ROM
This section includes the development of the 1D ROM for radial conduction in PCMs. It first includes the domain simplifications from the entire HX level to the segment-based level and assumptions. It then studies the condition required for the PCM to convert into a radial ROM. Lastly, it discusses the governing equations to simulate the ROM.

Domain Simplifications and Assumptions
Since the modeling of the entire HX domain is computationally expensive, it can be simplified to an independent segment-level-based model, using the following assumptions:

1.
Constant and isotropic material properties of fluid, polymer, metal, and PCM (both liquid and solid phases); 2.
Fully developed flow; Entire PCM at its melting temperature (T PCM ), initially; 5.
Heat transfer only through wire fins, and negligible conduction through polymer walls; 7.
Uniform mass flow rate ( . m total ) and temperature (T in ) profiles, per unit area at the inlet; 8.
Quasi-steady-state approximation (Ste < 0.1) or negligible thermal capacitance of PCM. Here, Ste = C p,PCM (T in − T PCM )/H. Since the fluid is assumed to have uniform mass-flux and temperature distribution across n f luid inlet channels in the x=0 plane of the entire TES domain (see Figure 2a), the fluid flow across every single channel would have the same mass-flow rate ( . m total /n f luid ) and inlet temperature (T in ). Now, using symmetry conditions, a single channel can be further reduced to a half-channel domain such that the wires of lengths W f luid /2 and W PCM /2 act as fins for the fluid and the PCM side, respectively. However, the single-row domain includes different fluid inlet temperatures ( for i = 1: ) that need to be calculated at a known time. As the fluid flows from any th wire to the 1 th wire, the fluid outlet temperature from th wire ( ) becomes the fluid inlet temperature ( ) for the 1 th wire and so on. In this way, for any th index of can be computed from the values of the previous wires. Since the This domain simplification results in an equivalent uniform mass-flow rate and uniform temperature of T in distributed over the entire half-channel in the y-direction. Now, for determining the inlet flow and temperature distribution across the wires in the tube-bank geometry, an index system of i (1 ≤ i ≤ N L ) in the x-direction and j (1 ≤ j ≤ N T ) in the y-direction was used to identify each wire. In this system, T i,j and . m i,j represent the flow temperature and mass-flow rate per unit area at upstream of the wire i, j, respectively and are called inlet temperature and inlet mass flow rate. As the temperature profiles are uniformly distributed across the inlet plane-x=0, the flow and temperature profiles for each ith wire would also be the same in the y-direction, i.e., T 1,j = T 1 for all 1 ≤ j ≤ N T (see Figure 2a). Extending this for any ith wire, T i,j = T i = f (y). Similarly, for the flow profile, as the mass-flow rate is uniformly distributed over the entire area across for any ith and jth wire, m. This means that for any given ith row, all the wires along the y-direction will have identical thermal performance. Thus, only N L wires need to be simulated instead of N T × N L . This simplifies the half-channel domain to a single-row domain, as shown in Figure 2a.
However, the single-row domain includes different N L fluid inlet temperatures (T i for i = 1:N L ) that need to be calculated at a known time. As the fluid flows from any ith wire to the (i + 1)th wire, the fluid outlet temperature from ith wire (T out i ) becomes the fluid inlet temperature (T i+1 ) for the (i + 1)th wire and so on. In this way, T i for any ith index of N L can be computed from the values of the previous wires. Since the fluid flow model across different N L wires is otherwise identical and independent, for a given T i , the model can be reduced from the single-row domain to a single PCM-wire Cartesian domain, as shown in Figure 2b.
However, the single PCM-wire Cartesian domain in Figure 2b includes three (Cartesian: x, y, z) directions of conduction for the PCM. As 1D radial conduction is assumed, the directions of conduction can be reduced from three to one (Radial: r). This simplification can be done by treating the single wire surrounded by the Cartesian PCM domain as a single wire surrounded by a cylindrical PCM domain. It contains the same thermal mass of the PCM as the Cartesian domain.
However, geometrically speaking, when the single-wire radial domains are arranged in the tube-bank geometry for all i and j values, some overlap regions might occur between the radial PCM domains of the neighboring wires (see Figure 3). However, geometrically speaking, when the single-wire radial domains are arranged in the tube-bank geometry for all and values, some overlap regions might occur between the radial PCM domains of the neighboring wires (see Figure 3). Thus, the next section entails the detailed analysis needed to compute and understand the geometric behavior of the overlap regions.

Domain Geometry of PCM-Wire Cylinders ( )
This section includes the realistic and simplified geometry consisting of the neighboring single-wire radial domains with the same thermal mass as the equivalent singlewire Cartesian domains (see Figure 4).
First, the maximum cylindrical radius is computed. Since the thermal mass is directly proportional to the PCM domain area for the same domain length-/2, the 2D PCM domain areas are computed for both the Cartesian ( ) and radial ( ) domains (see Figure 4) and are equated with each other to yield the condition for the maximum equivalent cylindrical radius ( ): where /2 and is the L.H.S. term while is the R.H.S. term. Thus, the next section entails the detailed analysis needed to compute r max and understand the geometric behavior of the overlap regions.

Domain Geometry of PCM-Wire Cylinders (r max )
This section includes the realistic and simplified geometry consisting of the neighboring single-wire radial domains with the same thermal mass as the equivalent single-wire Cartesian domains (see Figure 4).

Solving for
gives: (1) Thus, the maximum allowable PCM radius for the 1D ROM is obtained from Equation (1).  Now, as discussed, the PCM overlap regions and (see Figure 4b) formed due to the geometrical spacing of the adjacent cylinders, were studied. The overlap regions of the ROM domain could be problematic, as they can impact the conduction inside the region and the time taken to melt the PCM completely. Thus, these overlap areas between the cylinders, when represented in the Cartesian domain (see Figure 4a), are termed as gray areas. Here, the direction of conduction is affected by both of the neighboring PCM domains, which leads to 2D conduction in both x-and y-directions inside the PCM. In contrast, the equivalent overlap regions in the radial domains are assumed to have 1D radial conduction. Additionally, the time taken to completely melt the PCM in the radial domain would be relatively less than the time taken to completely melt the PCM in the Cartesian domain. This is because when the PCM in the overlap region is compared to its corresponding region in the Cartesian domain, the overlap region of the radial domain seems much closer to the wire than that of the Cartesian domain. This provides a shorter conduction path for the radial domain than the Cartesian domain and thus less melting time.
Mathematically, the overlap regions exist in five different cases controlled by the geometrical spacings between the neighboring PCM cylinders: (i) in the direction of (ydirection), (ii) in the direction of (see Figures 4a and 5 on xy plane), and (iii) in the direction of (x-direction). These different cases of overlap regions are governed by specified ranges of / ratios, which are derived mathematically (see Appendix A for complete derivation and analysis) and summarized as shown in Figure 6. Here, the controlled directions for the different overlap cases are also marked. It is to be noted that the neighboring PCM cylinders only cause these overlaps without overlapping the wires and thus are investigated only for that particular range of / . First, the maximum cylindrical radius is computed. Since the thermal mass is directly proportional to the PCM domain area for the same domain length-W PCM /2, the 2D PCM domain areas are computed for both the Cartesian (A cartesian ) and radial (A radial ) domains (see Figure 4) and are equated with each other to yield the condition for the maximum equivalent cylindrical radius (r max ): Solving for r max gives: Thus, the maximum allowable PCM radius for the 1D ROM is obtained from Equation (1). Now, as discussed, the PCM overlap regions A 1 and A 2 (see Figure 4b) formed due to the geometrical spacing of the adjacent cylinders, were studied. The overlap regions of the ROM domain could be problematic, as they can impact the conduction inside the region and the time taken to melt the PCM completely. Thus, these overlap areas between the cylinders, when represented in the Cartesian domain (see Figure 4a), are termed as gray areas. Here, the direction of conduction is affected by both of the neighboring PCM domains, which leads to 2D conduction in both x-and y-directions inside the PCM. In contrast, the equivalent overlap regions in the radial domains are assumed to have 1D radial conduction. Additionally, the time taken to completely melt the PCM in the radial domain would be relatively less than the time taken to completely melt the PCM in the Cartesian domain. This is because when the PCM in the overlap region is compared to its corresponding region in the Cartesian domain, the overlap region of the radial domain seems much closer to the wire than that of the Cartesian domain. This provides a shorter conduction path for the radial domain than the Cartesian domain and thus less melting time.
Mathematically, the overlap regions exist in five different cases controlled by the geometrical spacings between the neighboring PCM cylinders: (i) in the direction of S T (ydirection), (ii) in the direction of S D (see Figures 4a and 5 on xy plane), and (iii) in the direction of S L (x-direction). These different cases of overlap regions are governed by specified ranges of S T /S L ratios, which are derived mathematically (see Appendix A for complete derivation and analysis) and summarized as shown in Figure 6. Here, the controlled directions for the different overlap cases are also marked. It is to be noted that the neighboring PCM cylinders only cause these overlaps without overlapping the wires and thus are investigated only for that particular range of S T /S L .   The geometric behavior of the overlap regions is quantified by the ratio (A ratio ) of the PCM overlap region (A overlap = 2A 1 + A 2 ) to the total PCM region (A radial ). A detailed mathematical analysis is carried out in Appendix B. Thus, the geometric parameter A ratio can be determined by knowing the S T /S L and r max /r 0 values, as obtained from Equation (A12) in Appendix B and as shown below.

Governing Equations
As the TES design is set up, the governing equations of an analytical model of the TES are discussed. First, the modeling starts with the single-wire radial model to compute the phase-front location (r i (t)) for a given time-step by utilizing the thermal resistance approach based on energy conservation. In this way, the latent energy stored in an ith single-wire model (∆e i (t)) is determined. Second, the latent energy stored in the single-wire model is extended to the single-row domain (see Figure 2). Since the inlet temperatures are different for N L wire-rows (see Figure 2), each wire-row would have a different magnitude of the latent energy stored. Thus, all the wire-rows are individually iterated using the thermal resistance analogy, and the total latent energy stored (∆e total (t)) in the single-row domain (see Figure 2) is additively obtained. Third, the single-row model is extended to the entire TES model. As the single-row domain is identical for different PCM channels of the entire TES, the total latent energy, ∆E(t), stored in the entire TES is again additively obtained. Lastly, the total latent energy stored in the entire TES domain is integrated for different time-steps by predicting the next phase-front location (r i (t)) for the next time-step (i + 1). The model is simulated until the entire PCM has been completely transitioned or when latent energy is required for a specified period. Since each wire-row has a different magnitude of latent energy stored, the phase-front advancement will be different for each wire-row at different times. The wires closer to the inlet plane of the TES would completely melt more quickly than the ones placed farther away. The entire model was programmed in MATLAB.
The literature's traditional models for phase-change processes typically involve nonlinear transient heat equations [17], which are often expressed in terms of Stefan's number [37]. In general, Ste is assumed to be very small (less than 0.1), which means that the sensible energy is negligible as compared to the latent energy and thus a quasi-steady-state approximation is considered valid [17]. As the entire PCM is initially at its melting temperature (T PCM ), there is zero conduction present through the solid PCM. Thus, the steady-state thermal modeling of the ith single PCM-wire radial domain (see Figure 2b) includes only the conduction through the transitioned or liquid PCM, which can be computed analytically.
The total resistance for the ith single PCM-wire radial domain, R i total , at a time t is computed using the thermal circuit diagram from Figure 2b, as shown below. where the index i takes the value of 1 to N L wires. R i total can be obtained from the three resistances, which are connected in series across a known inlet temperature, T i , and PCM melting temperature, T PCM . These resistances are individually calculated for an ith single PCM-wire, as discussed in the following sections.

Calculation of Fluid Convective Resistance, R f luid
The fluid thermal resistance, R f luid , is calculated using Equation (4). where Here, h f luid is computed using CFD by considering flow parameters including Re number, flow properties including Pr number, and geometrical tube-bank parameters including S T D and S L D . Similarly, η f luid is calculated using h f luid , wire properties including k wire and geometrical parameters including W f luid and D. The authors previously studied and discussed this methodology as detailed in references [11,14]. Since the flow is assumed to be fully developed, h f luid is constant for all ith wires. This leads to constant η f luid and R f luid for all the ith wires (1 ≤ i ≤ N L ) and thus a steady-state condition is assumed.

Calculation of Wire Conductive Resistance, R cond
The conductive wire resistance of the portion of wire embedded in the polymer wall of thickness, t wall , (see Figure 2b) is calculated analytically [38] using Equation (5), by assuming 1D axial conduction.
where the cross-sectional wire area is A C = πD 2 /4. Since it is computed analytically and is independent of any ith wire row, it is considered a constant.

Calculation of PCM Conductive Resistance, R PCM (t)
Since the single PCM-wire radial model (see Figure 4b) assumes radial conduction inside the PCM, it treats the overlap region similarly to the non-overlap region and considers the modeling domain as a single entity of the PCM cylinder. The instantaneous PCM conductive resistance, R PCM , for any ith wire with a known melt-front location r i at time t is computed analytically using Equation (6).
where the index i takes the value of 1 to N L wires and A sur f PCM = πD W PCM 2 . This approach saves computational time compared to the otherwise transient and multi-dimensional CFD modeling of the PCM domain, as R i PCM is computed analytically. Here, the fin efficiency η i PCM for the single PCM-wire-domain is computed using a similar methodology used for the fluid side. The η i PCM for the wire-length of (L c PCM = W PCM 2 ) on the PCM side is obtained using the following Equations (7) and (8): where the index i takes the value of 1 to N L wires. Here, the fin parameter, b PCM , is defined as: where the index i takes the value of 1 to N L wires.
Since quasi-steady-state approximation is valid, h i PCM is computed analytically by assuming steady-state heat transfer through a cylindrical pipe [38]. It uses geometric parameters including r i and r 0 and PCM properties including k PCM , as shown in Equations (9) and (10).
where the index i takes the value of 1 to N L wires. Using Equations (4)- (6), R i total (t) is obtained (see Equation (3)) and is now used to compute the nominal heat transfer rate (q i 0 (t)), which can give latent energy (∆e i ) for the ith single PCM-wire at time t. The nominal heat transfer rate is computed using Equation (11).
where the index i takes the value of 1 to N L wires. However, in practice, we need to limit q i 0 (t) so that it is not larger than the maximum available heat transfer rate required to completely melt the available non-transitioned PCM (∆m max ) (see Figure 2b). Thus, the maximum available heat transfer rate (q max ) obtained from the maximum available latent heat related to ∆m max in ∆t time is evaluated as shown below.
where the index i takes the value of 1 to N L wires.
Here, ∆m i max (t) is computed based on the present melt-front location (r i (t)) and the maximum melt-front location (r max ), as shown below: where the index i takes the value of 1 to N L wires. Thus, the heat transfer rate of the ith wire can be calculated as shown in Equation (14), which ensures that the heat transfer for ith PCM-wire will be terminated once all surrounding PCM is melted.
where the index i takes the value of 1 to N L wires. Eventually, the latent energy stored by the ith single PCM-wire during this time-step is given in Equation (15).
where the index i takes the value of 1 to N L wires.

Latent Energy Stored in a Single Row, ∆e total (t)
The latent energy stored in the ith single PCM-wire is now integrated to compute the latent energy stored in an entire single row (see Figure 2b), consisting of a different total Energies 2022, 15, 5124 12 of 30 N L wires. The modeling of each ith wire can be determined using Section 3.3.1. However, the inlet temperatures for all the N L wires must be known (see Equation (11)). This can be obtained by using the energy balance between the nominal heat transfer rate for the PCM (q i (t)) and the sensible heat lost by the fluid across the ith wire, as shown below: where the index i takes the value of 1 to N L wires and . m is obtained from . m . As T i for all N L wires is now known, the latent energy stored (∆e total (t)) for all the N L wires is obtained using Equations (11)-(15) of Section 3.3.1. Thus, the total latent energy stored for the entire single row is obtained by adding the latent energy stored (∆e i (t)) for all the N L wires.

Latent Energy Stored in Entire TES, ∆E(t)
Similarly, the latent energy stored for the entire TES (∆E(t)) for the given time t is computed by simply summing up the latent energy stored in the single rows for the remaining N T rows of wires across n PCM channels of the HX.

Time-Integrated Latent Energy Stored in Entire TES, E total
As the latent energy stored in the TES is computed for a known r i for a given ith wire at time (t), the next step is to compute the latent energy stored in the TES for the subsequent time-steps. However, as the melt front progresses with time, an updated melt-front location, r i (t + ∆t), for a given ith wire, needs to be determined. Thus, the ROM uses the "forward Euler" method in time to compute r i (t + ∆t) for a given ith wire for the next discretized time-step, ∆t. As the latent energy storage for the single PCM-wire (∆e i (t)) leads to the melting of ∆m i (t) PCM (see Figure 7), the energy conservation at the melting interface is used to compute the updated melt-front location (r i (t + ∆t) for the next time-step, as shown in Equation (19).
where the index i takes the value of 1 to N L wires. As r i (t + ∆t) is known for each time-step, the time-integrated latent energy for the entire TES under a specified nominal time of t 0 can be obtained, as shown in Equation (21).
Graphically, typical trends obtained from the ROM for some given geometrical and operating conditions are shown in Figure 8. From Figure 8a, it can be seen that the heat transfer rate follows two trends. For t < 4.2 s, the heat transfer rate decreases continuously with time as the PCM resistance increases with time (see Equation (11)). Physically, as more PCM is melted with time, the heat travels a larger distance from the wire to reach the melt-front location, thus leading to higher PCM resistance and decreasing the heat transfer rate with time. next time-step, as shown in Equation (19).
where the index takes the value of 1 to wires. As Δ is known for each time-step, the time-integrated latent energy for the entire TES under a specified nominal time of can be obtained, as shown in Equation (21).
Graphically, typical trends obtained from the ROM for some given geometrical and operating conditions are shown in Figure 8. From Figure 8a, it can be seen that the heat   (11)). Physically as more PCM is melted with time, the heat travels a larger distance from the wire to reach the melt-front location, thus leading to higher PCM resistance and decreasing the hea transfer rate with time.
However, for 4.2 s, the heat transfer rate decreases discretely with time (see Equation (12)), where every stepwise drop corresponds to the complete melting of PCM of the th wire in the ascending order. The order of melting of the wires is better ex plained in Figure 8b, where it can be seen that there is more melted PCM in the first few rows ( 1,2,3) as compared to the last wire ( ), which is due to the decreasing ap proach temperature ( ) from 1 to . Furthermore, the area under the curve represents the total latent energy stored in the TES, . Thus, from a design perspective for a given stored latent energy, a TES that has a high average heat transfer rate and shor nominal time is advantageous.  In addition to the thermal performance of the TES, the hydrodynamic performance of the TES was computed using fluid side pressure drop (Δ ), determined via CFD However, for t > 4.2 s, the heat transfer rate decreases discretely with time (see Equation (12)), where every stepwise drop corresponds to the complete melting of PCM of the ith wire in the ascending order. The order of melting of the wires is better explained in Figure 8b, where it can be seen that there is more melted PCM in the first few rows (i = 1, 2, 3) as compared to the last wire (i = N L ), which is due to the decreasing approach temperature (T i − T PCM ) from i = 1 to N L . Furthermore, the area under the curve represents the total latent energy stored in the TES, E total . Thus, from a design perspective for a given stored latent energy, a TES that has a high average heat transfer rate and short nominal time is advantageous. In addition to the thermal performance of the TES, the hydrodynamic performance of the TES was computed using fluid side pressure drop (∆P f luid ), determined via CFD using flow parameters including . m and geometric parameters including S T D , n f luid , N L and S L D . The details of this process can be found in reference [11].

Non-Dimensionalized Form of Governing Equations
The governing equations of the 1D ROM are non-dimensionalized to obtain critical parameters dictating the TES performance. Here, the boundary condition is the major convective boundary on the wire exposed to the fluid side, with constant heat transfer coefficient h f luid (see Figures 2b and 9) and fluid-inlet temperature T i , while on the PCM side, the entire PCM that is not transitioned is assumed to be constant at its melting temperature, T PCM , continuously. Now, the fluid convective boundary condition in the 1D ROM can be further simplified such that the simplified fluid heat transfer coefficient h f luid yields the same resistance (1/(h f luid πr 2 0 )) as mentioned in Section 3.3 (R f luid + R cond ) for a given PCM mass. . The details of this process can be found in reference [11].

Non-Dimensionalized Form of Governing Equations
The governing equations of the 1D ROM are non-dimensionalized to obtain critical parameters dictating the TES performance. Here, the boundary condition is the major convective boundary on the wire exposed to the fluid side, with constant heat transfer coefficient h (see Figures 2b and 9) and fluid-inlet temperature , while on the PCM side, the entire PCM that is not transitioned is assumed to be constant at its melting temperature, TPCM, continuously. Now, the fluid convective boundary condition in the 1D ROM can be further simplified such that the simplified fluid heat transfer coefficient ℎ′ yields the same resistance (1/ ℎ )) as mentioned in Section 3.3 ( ) for a given PCM mass. The non-dimensionalization is performed by non-dimensionalizing the independent variables of the main governing equations, and into * and * , respectively, as shown in Equations (21) Here, , • , and * are the input parameters for the 1D ROM, as defined in Table 1. Table 1. List of non-dimensional groups.

Non-Dimensional Groups Mathematical Definition Physical Meaning
• Time-constant Proportional to the ratio of axial conductive resistance of the wire embedded in PCM to fluid-side resistance * 2 Proportional to the ratio of axial conductive resistance of the wire embedded in PCM to maximum radial resistance of PCM The non-dimensionalization is performed by non-dimensionalizing the independent variables of the main governing equations, r i and t into r * i and t * , respectively, as shown in Equations (21) and (22): Thus, the discretized governing equations of the 1D ROM (see Equation (A13) in Appendix C for details) can be non-dimensionalized as r i * and t* as derived in Equation (A21) in Appendix C and as shown below.
Here, τ, Bi·LR, and R * wire are the input parameters for the 1D ROM, as defined in Table 1.

Model Input and Outputs
The inputs for the model typically include the geometrical parameters, operating parameters, and thermophysical properties. The geometrical parameters include the TES design variables, such as wire diameter (D), wire spacings (S T and S L ), width of the PCM (W PCM ), and fluid channels (W f luid ). The operating parameters include fluid flow rate ( . m) and fluid-inlet temperature (T in ). Additionally, the thermophysical properties include the properties of air, metal wires, polymer, and PCM. These model inputs can also be grouped together into non-dimensionalized parameters as shown in Table 1.

Non-Dimensional Groups
Mathematical Definition Physical Meaning Proportional to the ratio of axial conductive resistance of the wire embedded in PCM to fluid-side resistance Proportional to the ratio of axial conductive resistance of the wire embedded in PCM to maximum radial resistance of PCM The outputs of the model include performance parameters, such as the melt-front location (r i (t)), heat transfer rate (q i (t)), total energy stored (E total ), and time taken to melt a specific quantity of the PCM.

2D CFD Cartesian Model (r, θ Plane)
This model was used to study the effect of the overlap areas between the neighboring cylinders (see Figure 4) on the performance parameters of the 1D ROM. Since the 1D ROM assumes 1D radial conduction, it treats the overlap region similarly to the non-overlap region. Thus, the model avoids the actual modeling of the 2D conduction of PCM present in the gray areas of the cartesian domain. Thus, the 2D CFD Cartesian model, which models the gray areas (see Figure 4b), was used to compare its performance parameters with that of the 1D ROM for a wide range of S T /S L and r * max . The error between the models is because the overlap region treated by the 1D ROM doesn't represent those gray areas of the Cartesian domain accurately. The present study considers the output performance parameter as non-dimensionalized time, t * , when 90% of the PCM has completely melted or when the liquid-fraction (χ) is 0.9. However, the discussion stands for other values of liquid fraction. Before simulating the 1D ROM and 2D CFD model, a condition of R * wire = 0 was considered for simplicity. This condition means that the wire has infinite thermal conductivity, or η PCM = 1, which means that there is uniform temperature along the entire length of the wire embedded in the PCM along the z-direction (see Figure 2b). This case also leads to uniform melting of the PCM in the z-direction. Thus, the cylindrical domain (r, θ, z) can be reduced into a 2D radial domain (r, θ plane) (see Figure 2b).
The radial domain (see Figure 4b) was simulated using the 1D ROM. Now, using the given condition on R * wire , the governing equation of the 1D ROM algorithm (see Equation (23)) yields a closed-form solution obtained in Appendix D and as shown below.
where R * f luid = R * wire Bi·LR (see Appendix D). The Cartesian domain (see Figure 4a) was simulated using a commercial CFD-based solver, Ansys Fluent, using standard enthalpy-based methods as detailed in the literature [39]. The CFD model comprises similar boundary and initial conditions as the 1D ROM (see Section 3.3.1). Here, a convective boundary of h f is imposed on the wires of cross-sectional area (πr 2 0 ), such that its convective resistance (1/ h f πr 2 0 ) is equivalent to the sum of convective resistances of the fluid side and wire (R f luid + R cond ) acting on the PCM side of the 1D ROM for the same PCM mass. The rest of the sides were taken as an adiabatic wall-boundary condition. The initial condition for the PCM was taken to be at its melting point (T PCM ). The rest of the fluid, wire, and PCM properties were taken to be similar to those in the 1D ROM.
An illustration of the trend of r * vs. t * for one such geometrical range of S T /S L =0.89 and r * max = 1.5 for conditions τ = 0.409 and R * f luid = 0.0518 is shown in Figure 10. Here, t * cartesian > t * radial , which means that, in the existing Cartesian domain, it takes a longer time for the PCM to be melted as compared to the simplified radial domain due to a longer conduction path. This could eventually affect the estimation of the latent energy of the TES at χ = 0.9. Thus, the overlapped effect is worth further study.

2D CFD Axisymmetric Model ( , Plane)
A 2D CFD axisymmetric model, typically used in the literature for cylindrical geometries [22,28], was used as a reference model to validate the performance parameters of the 1D ROM. The CFD model was simulated in Ansys Fluent, which utilizes an enthalpybased model [39]. The domain and boundary conditions of the 2D axisymmetric model are shown in Figure 11. Similarly, here, the convective boundary of ℎ was applied, as discussed in Section 3.4. and an adiabatic wall-boundary condition was applied to the rest of the sides. The 2D CFD axisymmetric reference model uses the same initial (see Section 3.3.1) boundary and operating conditions as the 1D ROM. However, the reference model includes both radial and axial conduction inside the PCM, while the 1D ROM assumes only radial conduction inside the PCM (see Section 3.1). Figure 11. 2D axisymmetric reference model.

2D CFD Axisymmetric Model (r, z Plane)
A 2D CFD axisymmetric model, typically used in the literature for cylindrical geometries [22,28], was used as a reference model to validate the performance parameters of the 1D ROM. The CFD model was simulated in Ansys Fluent, which utilizes an enthalpy-based model [39]. The domain and boundary conditions of the 2D axisymmetric model are shown in Figure 11. Similarly, here, the convective boundary of h f was applied, as discussed in Section 3.4. and an adiabatic wall-boundary condition was applied to the rest of the sides. The 2D CFD axisymmetric reference model uses the same initial (see Section 3.3.1) boundary and operating conditions as the 1D ROM. However, the reference model includes both radial and axial conduction inside the PCM, while the 1D ROM assumes only radial conduction inside the PCM (see Section 3.1).

2D CFD Axisymmetric Model ( , Plane)
A 2D CFD axisymmetric model, typically used in the literature for cylindrical geometries [22,28], was used as a reference model to validate the performance parameters of the 1D ROM. The CFD model was simulated in Ansys Fluent, which utilizes an enthalpybased model [39]. The domain and boundary conditions of the 2D axisymmetric model are shown in Figure 11. Similarly, here, the convective boundary of ℎ was applied, as discussed in Section 3.4. and an adiabatic wall-boundary condition was applied to the rest of the sides. The 2D CFD axisymmetric reference model uses the same initial (see Section 3.3.1) boundary and operating conditions as the 1D ROM. However, the reference model includes both radial and axial conduction inside the PCM, while the 1D ROM assumes only radial conduction inside the PCM (see Section 3.1). Figure 11. 2D axisymmetric reference model.

Test-Case Setup
The test case was set up by selecting an appropriate domain geometry, first studying Figure 11. 2D axisymmetric reference model.

Test-Case Setup
The test case was set up by selecting an appropriate domain geometry, first studying the geometric behavior of overlaps and then studying their effect on the performance parameters of 1D ROM.
A graph showing a variation of overlapped area ratio (A ratio ) versus the transverse to the longitudinal wire pitch ratio (S T /S L ), obtained from Equation (2), is plotted in Figure 12a. From the graph, it can be seen that for a given r * max , A ratio, or the overlap effect, varies across five different regimes bifurcated by the S T /S L values, corresponding to those shown in Figure 6. In the first regime (0.32 < S T /S L < 0.97), the overlap effect (A ratio ) starts decreasing with increasing S T /S L values as the neighboring PCM cylinders start moving further from each other, which decreases the common overlap region between them in the direction of S T . In the second regime (0.97 < S T /S L < 1.27), the overlap effect (A ratio ) continues decreasing in the direction of S T , but the overlap region in the direction of S D starts increasing as the PCM cylinders in the direction of S D start getting closer to each other. Thus, there exists a minimum overlap region of 1.27 around S T /S L . As S T /S L increases further, it transitions into the third regime (1.27 < S T /S L < 3.14) where the PCM cylinders in the direction of S T are much further from each other, yielding no overlap in the direction of S T . However, this brings the PCM cylinders closer to each other in the direction of S D , increasing the overlap region or A ratio . However, as S T /S L continues increasing in the fourth regime (3.142 < S T /S L < 4.12), the PCM cylinders now start distancing away from each other in the direction of S D , leading to a decrease in A ratio . However, the cylinders in the direction of S L start getting closer to each other, leading to an increase in A ratio . Thus, there again exists a minimum overlap region close to 3.5 around S T /S L . Finally, in the last regime, as the S T /S L increases beyond 4.12, there is no more overlap present due to cylinders in the direction of S D , but the cylinders in the direction of S L start getting closer to each other, resulting in increasing A ratio . Similarly, at a given S T /S L value, the A ratio is observed to decrease with increasing r * max because for a given overlap region (A overlap ), the available PCM domain area (A radial ) starts increasing. The minimum A ratio would eventually occur when r * max = ∞. As the 1D ROM doesn't capture the overlap effect, a detailed analysis was carried out to investigate the overlap effect on the performance parameters of the 1D ROM (t * radial (χ = 0.9)) by comparing it with that of the 2D CFD Cartesian model (t * cartesian (χ = 0.9)). The comparison between the models, %∆ = 1 − t * radial (χ=0.9) t * cartesian (χ=0.9) , was done for wide ranges of S T S L (0.32 to 12.57) and r * max (1.5 to ∞) (see Figure 12b). Graphically, it was observed that the deviation of the 1D ROM from the 2D Cartesian model varies across the five different regimes bifurcated by the similar S T S L values as in the case of A ratio . Thus, the %∆ vs. S T /S L trends seem to correlate well with the trend observed for A ratio vs. S T /S L . For S T /S L values with larger overlap effects (A ratio ), it would take longer melting times (t * cartesian ) in the case of the 2D Cartesian model due to a longer conduction path, than in the case of the 1D ROM (t * radial ). Longer melting times can lead to larger percentage deviation (%∆), i.e., t * cartesian > t * radial , as observed in Figure 10. The larger percentage deviation (%∆) or larger overlap ratios would be problematic for the 1D ROM in terms of TES design because it indicates the existence of high thermal resistance regions in the PCM and should thus be avoided. Similarly, smaller overlap regions (A ratio ) for the 2D Cartsian model would take almost comparable melting times as those predicted by the 1D ROM and, thus, smaller %∆ values. However, it should be noted that the calculated percentage deviation (%∆) is sensitive to mesh density and time-step.  From the above discussion, a test-case scenario was set up with S T /S L = 1.2 and r * max = 3, such that the overlap effect (A ratio ) and percentage error (%∆) are the mini- mum. The summarized geometrical and physical parameters for the test-case scenario are specified in Table 2. Since the 1D ROM utilizes a time-discretized approach, a time-step independence study was carried out to obtain time-step independent results. A typical time-step independence result for the test-case scenario involves quantifying the difference (%Diff) between the performance parameters, t * radial (χ = 0.9), between consecutive time-step sizes. The consecutive time-step sizes are obtained by refining the reference step-size by halving, as shown in Table 3. The performance parameter is reported to be time-step independent to within 1% when the step-size is halved, as the lowest time-step is sufficient enough to capture the simplest physics phenomenon of linear trend of r with respect to time. The equivalent number of steps (t 0 /∆t) corresponding to this reference are 300 and thus, this can be used as a reference for all the studies with a given t 0 .

Validation with the Widely-Used 2D CFD Axisymmetric Model
Since the 1D ROM is based on a segment-level single PCM-wire domain, the single PCM-wire model was validated with a reference 2D axisymmetric model (see Figure 11) based on test-case conditions (see Table 2). Here, the 2D axisymmetric model was simulated using a detailed parametric study as described in Section 5.

Parametric Case Study
The 1D radial model was validated with the reference 2D axisymmetric model for a wide range of dimensionless parameters, τ, Bi × LR, R * wire , as shown below, which act as the input parameters for the 1D ROM.
The output parameter for the current study is t * at χ = 0.9.

Effect of τ
The parameter τ, representing a dimensionless time-constant, is varied within the range specified in Table 4, while the other parameters are kept constant at the test-case conditions (see Table 2). As τ increases, which is also equivalent to the approach temperature increasing, there is more nominal heat provided to the PCM at a specific time, resulting in faster melting of the PCM or decreasing t * , as shown in Figure 13a. It can also be mathematically observed that τ is inversely related to t * (see Equation (23)). The ROM model also seems to be well-validated with the 2D reference model within 6%. Lastly, the validation of the 1D ROM with the 2D CFD axisymmetric model shows that the performance parameters * match well with each other within 3% for all cases * 1. However, for larger values of * , considerable deviations of about 15-65% were observed, which are more pronounced for smaller * values. This is mainly due to axial conduction occurring inside the PCM, which is not captured by the 1D ROM and is based on 1D radial conduction inside the PCM.

Effect of Bi·LR
The dimensionless parameter Bi·LR, proportional to resistance ratios of R wire to R f luid , is varied within the range specified in Table 4, while the other parameters are kept constant at the test-case conditions (see Table 2). An increase in Bi·LR for a given τ can also be expressed in terms of decreasing R f luid . As R f luid decreases, the total resistance also decreases for a constant R max PCM (see Equation (11)). This decrease in total resistance leads to an increase in the nominal heat transfer rate (q 0 ) provided to the PCM, resulting in faster melting of the PCM. Thus, t * decreases with increasing Bi·LR (see Figure 13b). This trend can be explained in three different cases.
In the first case, 0.03 ≤ Bi·LR ≤ 3. As Bi·LR increases from 0.03 to 3, t * suddenly decreases, mainly due to the sudden decrease in the total resistance. When Bi·LR = 0.03, the steady-state fluid-side resistance (R f luid ) is significantly higher and more-dominating than the transient PCM side (R f luid R max PCM ). This results in a significantly less nominal and almost steady-state heat transfer rate to the PCM, which takes a long time (t * ) to melt 90% of the PCM (see Equations (3) and (11)). In contrast, when Bi·LR increases from 0.03 to 3, R f luid decreases, but the transient PCM-side resistance now becomes slightly higher than R f luid , resulting in a sudden decrease in total resistance. Thus, a higher transient heat transfer rate occurs in the PCM and leads to a relatively shorter time (t * ) to melt 90% of the PCM.
In the second case, 3 ≤ Bi·LR ≤ 30. As Bi·LR continues increasing from 3 to 30, the fluid-side resistance continues decreasing for the same R max PCM such that R f luid /R max PCM <1. This results in a gradual increase in the transient heat transfer rate to the PCM, which leads to faster melting times (t * ).
In the last case, 30 ≤ Bi·LR ≤ 300. As Bi·LR continues increasing further from 30 to 300, the fluid-side resistance has already reduced so much that now, the equivalent total resistance doesn't change much and is nearly equivalent to the constant R max PCM when Bi·LR = 30. Thus, the melting times are almost the same for larger Bi·LR values for the 1D ROM.
Lastly, regarding the validation of the 1D ROM with the 2D CFD axisymmetric model, the performance parameters, t * , match well with each other within 7% for all cases Bi·LR ≤ 30. However, for Bi·LR around 300, more significant deviations up to 22% are observed, which are more pronounced for smaller values of t * . Though R f luid or R total decreases with increasing Bi·LR, for extremely large values of Bi·LR, the R wire becomes dominant in contrast to R max PCM , which results in significant temperature variations along the wire length and uneven radial melt-front advancement at different axial locations. Thus, the PCM starts conducting axially, while the 1D ROM is based only on the radial PCM conduction.
Effect of R * wire Similarly, the dimensionless parameter R * wire is varied within the range specified in Table 4. Here, R * wire is proportional to the resistance ratios of wire resistance for the wire embedded in PCM to the maximum PCM radial resistance. As R * wire is varied, the other parameters are kept constant at the test-case conditions (see Table 2). A decrease in R * wire for a given τ can also be expressed in terms of relative decrease in R wire or increase in k wire , which can be further expressed as an increase in fin efficiency on the PCM side (η PCM ). This decreasing R * wire leads to a decrease in PCM-side resistance R max PCM as well as a decrease in the total resistance, R total (see Equations (3), (6) and (11)). The decrease in R total results in an increase in the heat transfer rate provided to the PCM and the PCM starts melting faster. Thus, t * reduces with decreasing R * wire , as shown in Figure 13c, for the 1D ROM. This trend can be explained in three different cases.
In the first case, 100 ≥ R * wire ≥ 10. When R * wire = 100, the R wire becomes dominating in contrast to the R max PCM (R wire R max PCM ), resulting in non-uniform temperature distribution present along the wire length on the PCM side. This creates significant temperature variations along the wire length and uneven radial melt-front advancement at different axial locations, leading to more axial conduction inside the PCM than radial. Thus, PCM takes more time to completely melt radially. As R * wire decreases from 100 to 10, the R wire to R max PCM suddenly reduces by nearly six times. Thus, the axial conduction inside the PCM is less dominating, and there would be slightly more radial conduction inside the PCM, leading to reduced time to melt the PCM completely radially, t * .
In the second case, 10 ≥ R * wire ≥ 1. When R * wire = 1, the R wire becomes almost comparable to the fluid-side resistance, leading to much more radial conduction inside the PCM, and thus t * further reduces.
Similarly, in the third case, 1 ≥ R * wire ≥ 0.01. As R * wire further reduces from 1 to 0.01, the R wire continues to decrease, but the fin efficiency increases only slightly as it reaches its limiting value of 1. Thus, the equivalent R max PCM doesn't decrease much from when R * wire = 1. This leads to a minimal change in R total and thus a minimal change in the nominal heat transfer rate to the PCM, as per the 1D ROM. Thus, it takes almost the same time to melt the PCM completely as when R * wire = 1. Lastly, the validation of the 1D ROM with the 2D CFD axisymmetric model shows that the performance parameters t * match well with each other within 3% for all cases R * wire ≤ 1. However, for larger values of R * wire , considerable deviations of about 15-65% were observed, which are more pronounced for smaller t * values. This is mainly due to axial conduction occurring inside the PCM, which is not captured by the 1D ROM and is based on 1D radial conduction inside the PCM.

Multi-Dimensional Parametric Study
The model's validity over a wide range of dimensionless parameters was studied by simultaneously varying the dimensionless parameters for each possible combination, as seen in Figure 14. The 1D ROM seems to correlate well with the 2D CFD axisymmetric model within 10% for most of the values of Bi·LR and R * wire , except for a few extreme cases. When Bi·LR ≥ 300, for extremely high values, the R wire is significantly higher than R max PCM , resulting in only dominant axial conduction inside the PCM. Since the 1D ROM is based on 1D radial conduction inside the PCM, this effect isn't captured, so there is considerable deviation of about 86% in t * at χ = 0.9 between the 1D ROM and the 2D CFD axisymmetric model.

Multi-Dimensional Parametric Study
The model's validity over a wide range of dimensionless parameters was studied by simultaneously varying the dimensionless parameters for each possible combination, as seen in Figure 14. The 1D ROM seems to correlate well with the 2D CFD axisymmetric model within 10% for most of the values of • and * , except for a few extreme cases. When • 300, for extremely high values, the is significantly higher than , resulting in only dominant axial conduction inside the PCM. Since the 1D ROM is based on 1D radial conduction inside the PCM, this effect isn't captured, so there is considerable deviation of about 86% in * at 0.9 between the 1D ROM and the 2D CFD axisymmetric model. Similarly, when * 10, for extremely high values, the is significantly higher than and thus leads to axial conduction inside the PCM and more significant deviations in * between the 1D ROM and 2D CFD axisymmetric model. However, there might be cases when ≫ : i.e., when • 3 such that the fluid-side resistance dominates and leads to almost constant steady-state total resistance, which happens in some applications. In such cases, the transient behavior of the heat transfer inside the PCM doesn't matter much, and the validation between the 1D ROM and 2D CFD axisymmetric model would always be unconditionally true. Thus, the wide ranges facilitating radial conduction in the 1D ROM model were successfully identified.

Conclusions
An efficient validated 1D ROM was formulated, which is analytically based on the thermal resistance analogy and conservation of energy, and which accurately models the performance of the TES studied here. The ROM was applied in characterizing the performance of a cross-media-based, fluid-to-PCM, metal-polymer composite HX, which is a Figure 14.
Comparison between 1D ROM and the 2D CFD axisymmetric model for performance parameter.
Similarly, when R * ≥ 10, for extremely high values, the R wire is significantly higher than R max PCM and thus leads to axial conduction inside the PCM and more significant deviations in t * between the 1D ROM and 2D CFD axisymmetric model. However, there might be cases when R f luid R wire : i.e., when Bi·LR ≤ 3 such that the fluid-side resistance dominates and leads to almost constant steady-state total resistance, which happens in some applications. In such cases, the transient behavior of the heat transfer inside the PCM doesn't matter much, and the validation between the 1D ROM and 2D CFD axisymmetric model would always be unconditionally true.
Thus, the wide ranges facilitating radial conduction in the 1D ROM model were successfully identified.

Conclusions
An efficient validated 1D ROM was formulated, which is analytically based on the thermal resistance analogy and conservation of energy, and which accurately models the performance of the TES studied here. The ROM was applied in characterizing the performance of a cross-media-based, fluid-to-PCM, metal-polymer composite HX, which is a novel HX design concept providing an in-built macro-encapsulation for the PCM. The advantages of the developed system include higher performance, lighter weight, and a lower cost compared to conventional HXs in the same class of applications. The 1D ROM model simplifies the entire HX domain to a segment-level domain by treating the PCM as cylindrically wrapped around the single wire, based on geometrical spacings related to the wires. A mathematical analysis was carried out to study the effect of these geometrical spacings on the validity of the ROM model. Optimum spacing ratios of S T /S L = 1.27 and 3.5 at higher r * max were identified for accurate results. The 1D ROM model was then subsequently validated with a 2D axisymmetric model, extensively used in the literature for cylindrical PCM studies. The validation was performed for a wide range of dimensionless parameters: (i) τ (ranging from 0.03 to 300), based on time constant, (ii) Bi·LR (ranging from 0.03 to 300), based on axial conductive resistance of the wire embedded in PCM to fluid-side resistance, and (iii) R * wire (ranging from 0.01 to 100), based on the axial conductive resistance of the wire embedded in PCM to maximum radial PCM resistance. The 1D ROM was found to be accurate within 10% compared to the 2D axisymmetric model for almost all the ranges except for extreme ranges when Bi·LR ≥ 300 or when R * wire ≥ 10, such that there is significantly large resistance in the wire embedded in the PCM, leading to axial conduction in the PCM. However, for very low values of Bi·LR ≤ 3, the fluid-side resistance becomes dominant, leading to a gradual steady-state heat transfer process independent of the R * wire value. Thus, the model would unconditionally correlate well with the 2D CFD axisymmetric model for such cases. For applications with a known τ value, requiring shorter melting times of PCM, the TES design for larger Bi·LR (<300) and smaller R * wire conditions (≤1) is recommended. In contrast, for applications requiring longer melting times of PCM, the TES design for smaller Bi·LR (<3) and larger R * wire conditions (around 1) is recommended. All the applications should be designed near optimum spacing ratios, with higher r * max . Thus, the validated and robust 1D ROM can be used to predict latent thermal energy storage for such PCM-fin-based geometries and for diverse applications without using the computationally expensive CFD-based tools.

Appendix A. Range of Overlaps
The range of S T S L values determine the nature of overlaps. There can be major five cases of overlap types depending on their geometrical spacing between the neighboring PCM cylinders.

•
Overlap in direction of S T or S T controlled For this to occur, S T < 2r max (see Figure 5a,b). Substituting r max from Equation (1) and simplifying the above equation in terms of S T /S L , we get: However, for S T S L < 1.27, there can also be overlaps between n-PCM cylinders kept at nS T distance apart, when n > 1. For that to occur, the nS T < 2r max condition should be satisfied, yielding: S T S L < 1.27 n 2 (A2) However, our present study is limited to overlaps only due to neighboring cylinders (n = 1). Thus, the equivalent the S T /S L range can be defined as: 0.32 < S T S L < 1.27 (A3) where the lower limit of S T /S L occurs when n = 2 in Equation (A1).
• Overlap in direction of S D or S D controlled Similarly, for this to occur (see Figure 5b-d), we need to meet the following conditions but instead of S T , we have S D : where S D = (S T /2) 2 + (S L ) 2 . The expression of Equation (A4) in terms of S T and S L and further simplification of the inequality yields the actual bounds for overlap between the neighboring PCM cylinders in the direction of S D as: 0.97 < S T S L < 4.12 (A5) Moreover, it is found that any overlaps possible between n-PCM cylinders (n > 1) aren't geometrically possible. Thus, Equation (A5) gives the actual bounds.

•
Overlap in direction ofS L or S L controlled For this to occur (see Figure 5d,e), we need the minimum distance between neighboring PCM cylinders in the longitudinal direction (2S L ) to be less than the diameter of the PCM cylinders of maximum radii. Further simplification of the inequality gives the following condition: S T S L > 3.142 (A6) Moreover, to avoid any possible overlaps between n-PCM cylinders, the new limits for S T S L would be: 3.142 < S T S L < 3.142n 2 (A7) Thus, the equivalent range is: 3.142 < S T S L < 12.6 (A8) where the upper limit of S T /S L occurs when n = 2 in Equation (A7) Lastly, to ensure that the PCM cylinders don't overlap with the wires, the following condition must be satisfied: r 0 ≤ min(r max , p − r max ) where p = min(S T , 2S L , S D ), valid for all types of overlaps in any distance. Graphically, summarizing the above overlap-types gives five possible cases, as shown in Figure 5.
where τ = Ste × Fo is a dimensionless constant and is defined in Table 1.
Here as typically used in the literature [15].
Fourth, Equation (A17) is now rearranged by simplifying the η PCM term. The η PCM term is simplified from the main source Equation (7), as shown in the following equation: where r * max = r max /r 0 . Let resistance-related parameters be defined as R max cylinder = log(r * max ) 2πk PCM W PCM and R wire = W PCM k wire πr 2 0 . Thus, Equation (A18) can be expressed as: where R * wire = R wire R max cylinder log(r * max ) is some dimensionless constant (see Table 1), which is also found to be equivalent to 2k * LR 2 .
Fifth, Equation (A17) is now rearranged by simplifying the ((R f luid + R cond ) (k PCM πW PCM )) term. As the computational domain of the single PCM-wire model (see Figure 2b) can be further simplified into Figure 9, Equation (A17) gets simplified to: τdt * = r * k wire h f luid (2r 0 ) Let the terms be defined in known dimensionless parameters present in the literature [15], Bi = h f luid (2r 0 ) k wire , and let there be simple ratios such as thermal-conductivity k * = k PCM k wire and LR = W PCM r 0 . Together, these terms can be grouped as Bi·LR, which is also found to be equivalent to 2R wire /R f luid (see Table 1).
Thus, Equation (A17) can be expressed as shown below: Appendix D. Limiting Conditions: R * wire = 0 A limiting case might arise when R * wire = 0 as k * = 0 delivering η PCM = 1. In such a scenario, Equation (A21) is modified by substituting R * wire = R wire R max cylinder log(r * max ) for 2k * ·LR 2 and Bi·LR = 2R wire /R f luid , as shown in the following equation: where R * f luid = R f luid R max cylinder log(r * max ) is similar to the definition of R * wire . Thus, Equation (A22) can be easily integrated from (1, 0) to (r * i , t * ) as shown below: