Dynamic modelling of ice ‐ based thermal energy storage for cooling applications

The development of accurate dynamic models of thermal energy storage (TES) units is important for their effective operation within cooling systems. This paper presents a one ‐ dimensional discretised dynamic model of an ice ‐ based TES tank. Simplicity and portability are key attributes of the presented model as they enable its implementation in any programing language which would, in turn, facilitate the simulation and analysis of complex cooling systems. The model considers three main components: energy balance, definition of the specific heat curve, and calculation of the overall heat transfer coefficient. An advantage of the model is that it can be adapted to other types of TES units employing phase change materials. The modelling approach assumes equal flow and temperature distribution in the tank and considers two internal tubes only to represent the whole tank—significantly reducing the number of equations required and thus the computation time. Thermophysical properties of water during the phase change and of the heat transfer fluid are captured. The ice ‐ based TES tank model has been implemented in MATLAB/Simulink. A good agreement between simulation results and experimental data available in the literature has been achieved—providing confidence in the validity of the mathematical model.


| INTRODUCTION
Cooling systems are used to provide comfortable air conditions in buildings and to deliver refrigeration and cooling services to manufacturing processes [1,2].Typically, a cooling system comprises an electric chiller and a heat exchanger.A hydraulic network with pumps, pipes and an arrangement of control valves is used to circulate a heat transfer fluid (HTF) to meet cooling demands.The chiller uses electricity and chilled water supplied from cooling towers to reduce the temperature of the HTF.
The time mismatch between low electricity costs (for chiller operation) and peak cooling demands and the complexity of the hydraulic circuits used in cooling networks create important challenges.To address these issues, thermal energy storage (TES) units can be incorporated into cooling systems to act as a buffer between supply and demand and to provide flexibility.This enables the peak cooling demand to be shaved, electrical load to be shifted and electricity costs reduced.For instance, a TES tank may act as an energy sink during minimum cooling loading, when cost of electricity is reduced, and operate as a cooling source during peak demand, when cost of electricity is high.
Phase change materials (PCMs) are commonly used as storage media in TES units [2][3][4].A PCM is a substance that absorbs or releases a large amount of energy during its phase transition-which occurs at a nearly constant temperature.This energy is referred to as the latent heat value (LHV) [5].A high LHV increases the storage capacity of PCM-based TES units compared to other storage media relying on temperature changes only (i.e., relying on sensible heat) [6].Although there are many types of PCM available, ice is a preferred choice for cooling applications due to its high energy density, low cost and, particularly, its melting temperature [7].
- 1     Existing mathematical models of PCM-based TES tanks consider the internal structure of the storage tank, the type of PCM and the HTF [6,[8][9][10].Steady-state modelling remains a widely used approach.Although energy balance (based on differential equations) is employed to describe the thermal behaviour inside the tank, the finite difference method using algebraic equations is adopted to iteratively solve the system [10][11][12][13][14][15][16][17][18][19][20].Thus, changes exhibited during real-time operation are disregarded.For instance, the heat transfer coefficient, which is key to determine the dynamic behaviour of a heat transfer process, is assumed constant.Moreover, specific heat or enthalpy curves, which comprise essential thermophysical properties describing the thermal behaviour of the PCM, are not considered.For instance, [11,12] use an algebraic equation based on the LHV to describe the enthalpy of the PCM during the phase transition.In [13] the product of the LHV by the total mass of the PCM is used to calculate the released or absorbed energy.Although some references provide useful methods for the approximate calculation of heat transfer coefficients, the effects of the thermophysical properties of the HTF are not accounted [14,17,18].
Three-dimensional (3-D) finite element-based models accurately describe the thermal dynamics inside the storage tank [18][19][20][21][22][23][24].This includes representing the behaviour of the PCM as a source of latent heat and of the HTF in terms of liquid fraction, energy transferred and momentum.Simulation outputs are commonly used to design the internal structure of the storage tank so that its heat transfer efficiency is improved.However, 3-D models require not only high-performance computing, but also large simulation times.
Two-dimensional (2-D) and one-dimensional (1-D) models are less complex than their 3-D counterparts.In these models, lumped control volumes are defined, where only variations through the x direction for 1-D representations and x and y directions for 2-D representations are considered.Given that changes through the z-direction are neglected, the time for any given solution shortens.Reference [25] presents a 1-D model of a heat exchanger with PCM.This work is significant as it shows the use of the specific heat curve to define the LHV of the PCM; however, the calculation of the heat transfer coefficient is not explained.The 1-D model of a PCM-based TES unit is presented in [26], where a basic equation considering a constant enthalpy represents the phase transition of the PCM.The accuracy of the models in [25,26] is arguably low when simulation results are compared to experimental data.
A better agreement between simulation and experimental results using 1-D models is shown in [27,28].In [27], the model of a high-grade cold storage unit with PCM is presented.An enthalpy curve is defined by a mathematical equation.However, an iterative process is used to solve the discretised model which, in turn, prevents its use in the analysis of more complex cooling systems.A PCM-based TES tank is modelled as a 1-D representation in [28] and experimental results are presented.In this reference, the LHV of the PCM is defined by a specific heat curve and an empirical function is used to compute the convection heat transfer due to the HTF.
With regards to 2-D models, a comparison between simulation results with experimental data is included in [29,30].An active air-PCM heat exchanger unit is presented in [29].Although there is a good agreement between simulation and experimental results, the enthalpy curve of the PCM is not shown.Conversely, a 2-D model of a PCM-based TES unit is presented in [30], where detailed calculation of the heat transfer coefficient and a well-defined specific heat-enthalpy curve of the PCM are considered.This work is relevant as the analysis allows to compute temperature variations using energy balance through the x and y directions.
This paper presents a dynamic yet simple 1-D mathematical model of an ice-based TES tank for cooling applications.The model is defined by a set of nonlinear differential equations and uses energy balance to describe the thermal behaviour inside the tank.It considers the thermophysical properties for the temperature range of interest of the water/ice (i.e. the PCM), of the HTF, and of the tubings enclosing the HTF.The heat transfer phenomenon between the HTF and the PCM is modelled dynamically.Unlike ice-based tank models found in references [31,32], it uses a specific heat curve instead of a constant LHV to dynamically describe the phase change.The model also considers the internal structure of the tank, including the heat transfer area, hydraulic diameter, thickness and length of the tubes.Even when the model encapsulates relevant hydraulic and thermal dynamics, it achieves a reduced computation time compared to complex 3-D models and it is simpler than 2-D models.
The tank model was implemented in MATLAB/Simulink to analyse different operating conditions, including charging and discharging processes.To verify the validity of the model, simulation results are compared against experimental results reported in [33,34], where a specific ice-based tank configuration is presented but no mathematical model of the storage tank is included.
The novelty of the presented model stems from its simplicity and portability.By considering a representative internal structure of the ice-based TES tank, the description of its thermal behaviour is simplified.However, the model is sufficiently detailed to accurately describe heat transfer during charging and discharging processes under different mass flow rates of the HTF.Moreover, it captures the temperature dependence of the thermophysical properties of the HTF and of water/ice.In addition, the model can be implemented in any programing language with an engine solver for ordinary differential equations (ODEs).It also considers a fitting process with experimental results.This is encapsulated by a correction factor for the calculation of the overall heat transfer coefficient.
The main contributions, strengths, and advantages of this work are summarised as follows: � The model is adaptable to different internal geometries of a TES system where the shape of the HTF channel is welldefined.For instance, slab and cylindrical containers of PCM would form rectangular and annular channels, respectively.
� The temperature dependency of the thermophysical properties of the HTF and the PCM is included in the mathematical model.Particularly, the specific heat curve is defined using probability density functions, enabling further accuracy of the dynamic model-particularly during phase change.By modifying the coefficients of these functions, the specific heat curve can be adapted to a different type of PCM (and extended to heating applications).� The mathematical model has been provided as a pseudocode, which may be implemented in any software package with an ODE solver.This makes the model portable and suitable for system-level studies of cooling energy centres or advanced industrial processes.� The simplicity of the mathematical model leads to a significantly reduced computation time compared to more complex representations.� The developed model has been verified against experimental data from a commercial ice-based tank available in the literature.Simulation and experimental results agree on well despite the model's simplicity-providing confidence in the modelling approach.To the knowledge of the authors, no other work providing a sufficient level of modelling detail presents such a comparison between simulation and experimental results for an ice-based TES tank.

| SYSTEM UNDER STUDY
Figure 1 shows a typical TES configuration within a cooling system energy centre.It includes two ice-based TES tanks.A valve arrangement enables switching between different operating modes.During the charging process, the HTF is cooled down to temperatures below 0°C by the compressor chiller and is pumped to the tanks.To achieve this, valve 4 is opened and valves 1, 2 and 3 are closed to allow the cold HTF to pass through the tanks.Conversely, the chiller is turned off for a discharging process, where valves 2 and 4 are closed and valves 1 and 3 are opened (and the flow is bypassed through a pipe and a valve before the compressor input).If the ice is completely depleted following a discharging process, a third operating mode is also possible to continue supplying cooling services without using the ice stores.In this case, the chiller is turned on, valve 1 remains open, valves 3 and 4 are closed, and valve 2 is opened.The ice-based TES tank considered in this work is a Calmac ICEBANK 1098C model.This has a nominal storage capacity of 350 kWh [35].Figure 2(a) shows a schematic for the interior of the tank.The arrangement consists of 34 spiralled pairs of polyethylene tubes (i.e.34 horizontal levels � 2 tubes = 68 tubes in total) submerged in water.Each pair of tubes a and b is rolled in a spiral shape and placed in a horizontal plane throughout the height of the tank, as shown in Figure 3. HTF at a suitable temperature is pumped inside the tubes to solidify water into ice (for charging) or melt the ice built around the tubes into water (for discharging).
As shown in Figure 2(b) (reproduced from [31]), there are four vertical tube headers that connect each horizontal level to the other levels.An internal header is connected to an external one to supply the HTF and the remaining internal and external headers are coupled for the return path of the HTF (i.e. the HTF inlet and outlet headers, respectively).This arrangement enables the HTF inside the pair of tubes for any horizontal level to flow in a counter-current direction to, in turn, facilitate homogenous ice growth and melting throughout the tank.In addition, the split of the flow through the 68 tubes reduces the pressure drop due to the HTF circulating through the tank.Figure 3 also shows a top view of a single horizontal level for the tank's charging and discharging modes.For a charging process, ice builds in a cylindrical shape around each tube when the HTF cooled down to low temperatures (e.g.−4°C) circulates inside the tubes.During a discharging process, a warm HTF (e.g.12-15°C) flows through the tubes so that the cooling energy stored in the ice is transferred to the HTF by melting the ice.
Figure 4 provides further detail on the ice formation for adjacent tubes during the charging process.The specific dimensions of the cylindrical shape of water/ice around the tubes is defined by the radius of the ice, r ice .This assumption helps establishing the limits of the control volume for the HTF, tubes and water/ice.This volume is used to apply energy balance to determine the average temperature of the HTF and of water/ice.In the figure, the tube dimension is defined by its internal (r i ) and external (r o ) radii.

| MODELLING APPROACH
A dynamic model of the tank which describes the thermal behaviour of the HTF and of the water/ice was developed.The heat transfer process inside the tank considers forced convection between the HTF and the internal wall of a pair of tubes, conduction in a radial direction into the tubes and the volume of water/ice, and conduction between the control volumes of the tubes.The phase change of water was modelled using the relationship between enthalpy and specific heat.
Two important assumptions are made: an equal internal flow and an equal temperature distribution for all pair of rolled tubes throughout the tank's height exist.Since every pair of tubes has the same conditions, the energy balance applied for two tubes only is representative of the behaviour for all tubes and, thus, the analysis is restricted to a single pair of tubes.The tank is assumed to be well-insulated; therefore, energy losses due to the ambient temperature are ignored.

| Energy balance
To apply energy balance to the ice-based tank, it is necessary to establish a control volume for each tube.For this paper, the control volume is defined by a tube section which includes the HTF volume inside the tube, the water/ice around it, and the wall of the tube through which the heat is transferred.This is shown in Figure 5 through a lateral and frontal crosssectional view of a single tube.For simplicity, the energy balance is done for a single tube, but the modelling approach and the analysis consider a pair of tubes.
The total thermal energy transferred to a specific element (e.g. the temperature changes of a tube wall due to hot fluid passing through it and the environment temperature) is given by the sum of the heat contributions from different sources (e.g.steam produced at certain temperature, or hot surface heating by electrical resistances).Dynamically, this transit of energy is described by where the overdot notation stands for a variation of a given variable (i.e.rate of change with respect to time), E [ J] is the energy stored, Q in [ J/s] is the heat provided by an external source and Q tr [ J/s] is the heat given by the energy transferred to/from another element.The scope of this work is restricted to a 1-D analysis.In this case, an energy balance for the HTF describes how the mean temperature T f [°C] of the fluid varies with respect to the x direction along the tube only and how the convection heat transferred is affected by this change.To further reduce model complexity, viscous dissipation, radial fluid flow, axial heat conduction, external forces, and compressibility are ignored [30,36,37].The heat transferred to a specific volume of HTF is described by the variation of its internal energy [37], which is defined by where E f [ J] is the internal energy of the HTF, m f [kg] its mass and c p,f [ J/kg°C] its specific heat, which is a function of temperature, d denotes a differential value, and dT f [°C] is the differential of the mean temperature of the HTF's mass.The total heat provided by the external source is related to the difference between the input and output temperatures through the differential control volume (see Figure 5).This is defined as [37].
where T f,in is the input temperature and _ m f [kg/s] is the mass flow rate of the HTF.Assuming the mean temperature of the control volume as T f,in + dT f , the heat convection within the control volume is expressed as [37].
The heat transferred Q tr between water/ice and the HTF depends on their temperature difference.It is calculated using the overall heat transfer coefficient U [W/m 2 °C] as follows: where A tr [m 2 ] is the heat transfer area and T w [°C] is the mean temperature of the volume of water/ice.The radial direction of the heat rate is illustrated by the red and blue arrows shown in Figure 6.Heat flowing from the HTF to the water/ice (red arrows) implies a higher temperature of the HTF, which melts the ice.Conversely, a higher temperature of water/ice causes heat to flow towards the HTF (blue arrows), which increases the HTF temperature during ice formation.It is assumed that heat is transferred into the cylindrical volume of water/ice by conduction in the radial direction only.The heat transfer between the ice/water control volumes of the tubes (Q w,a and Q w,b ) is also illustrated by the white and black arrows.The heat transfer phenomena (conduction and convection) are shown at the bottom of Figure 6 through an electrical analogy [38].
Using (1)-( 5), and considering the mass of the HTF in terms of its density ρ f [kg/m 3 ] and volume V f [m 3 ], the rate of change of energy _ E f ;a of the HTF for tube a is calculated by where subscript 'a' is used to refer to tube a.
With regards to the rate of change of the thermal energy of water/ice _ E w;a for tube a, it is assumed that the water is completely stagnant.This implies that no natural convection or radiation are exhibited in the liquid phase of the water/ice.Thus, there are two heat sources for tube a: the conduction heat transfer given by the interaction with the HTF and the conduction heat transfer with the control volume of tube b (see Figure 6).The first heat source has the same magnitude as the heat transferred in ( 6), but with an opposite direction.The second heat source is defined by where A e [m 2 ] is the external surface of the control volume, subscript 'b' is used to refer to tube b, T w,a and T w,b [°C] are the temperatures of water/ice of the control volume of tubes a and b, and U w [W/(m 2 °C)] is the conduction heat transfer coefficient between the control volumes of water/ice of the tubes, defined by where k w [W/(m°C)] is the thermal conductivity of water, r o [m] is the external radius of a tube, and r ice [m] is the assumed radius of the control volume of water/ice (see Figure 4 for reference).Although the exact contact area of the control volumes is not well defined due to the circular shape of the control volume, it is assumed that heat is exchanged through the whole surface area A e .An additional heat source term, _ E l , is introduced to account for the water/ice that may not melt or solidify during the discharging and charging processes.Given that the tank is assumed well-insulated, the ambient temperature is not considered and the heat loss results only from the water/ice remaining at melting temperature.Thus, the heat transferred is defined in terms of the temperature difference between the HTF and water/ice (Q tr ), and of the control volumes of water/ice of each tube (Q w,a , Q w,b ).By using the thermal energy variation integrated over the cross-sectional volume as for the HTF in (2), _ E w;a is described by where c p,w [ J/(kg°C)] is the specific heat of water/ice, V w [m 3 ] its volume, and ρ w [kg/m 3 ] its density.The heat loss _ E l;a for tube a is defined by where U l [W/(m 2 °C)] is the conduction heat transfer coefficient.As the heat loss is caused by the portion of water/ice whose phase change is incomplete, the second temperature term is kept at 0°C.The right-hand side terms of ( 6) and ( 9) are, respectively, the variations in internal energy of the HTF and of water/ice for tube a.In equation ( 9), the specific heat (c p,w ) is a key thermophysical property for the phase change of water to ice and vice versa.The equations for the energy balance of tube b have similar terms as ( 6) and (9).Subscript a is simply replaced by subscript b and an opposite direction of heat transfer between the control volumes of water/ice is considered: Note: Considering all the factors that may affect heat loss in the TES tank represents a challenging task.Usually, a lumped factor is assumed to simplify this process.Although an accurate calculation of U l is beyond the scope of this work, a heuristic approach was adopted to find its value as 0.4084 W/(m 2 °C).To this end, simulation and experimental results were compared and the value of U l was adjusted iteratively until a good agreement between the sets of results was achieved.This process is not cumbersome as all parameters in (10) besides U l are welldefined.

| Overall heat transfer coefficient
The overall heat transfer coefficients for the charging and discharging processes are typically provided by manufacturers or obtained experimentally, but they can also be calculated analytically considering the mechanical conditions of the HTF inside the tubes.In [31], given the difficulties to distinguish a change from laminar to turbulent flow in a curved tube, the Reynolds number is assumed as constant and only laminar flow is considered.However, turbulent flow is discussed in this paper given the wide range of mass flow rates associated to the ice tank operation.As opposed to [31], the fluid mechanics conditions inside the tubes and convection heat transfer are here considered to dynamically calculate the overall heat transfer coefficient.
Let the mean velocity v f [m/s] of the fluid through the internal cross-sectional area of the tube A c [m 2 ] be defined by where r i [m] is the internal radius of the tube (see Figure 4).The Reynolds number Re is calculated as is the dynamic viscosity of the fluid.The value of Re determines whether the fluid through the tube is laminar (Re ≤ 2300) or turbulent (Re > 5000).For the sake of simplicity, the transition zone between laminar and turbulent fluids is not considered.For laminar flow, the Nusselt number Nu, a key element for calculating convection heat transfer, remains constant, with a value of 4.36 [37].On the other hand, if the fluid is turbulent, Nu is calculated as where Pr is the Prandtl number of the HTF, n = 0.4 for a heating process and n = 0.3 for a cooling process.
Using (16), the convection heat transfer coefficient U cv [W/m 2 °C] of the HTF is calculated by where k f [W/m°C] is the thermal conductivity of the fluid.The conduction heat transfer coefficient in the tube wall U cd,t [W/m 2 °C] is calculated by where k t [W/m°C] is the thermal conductivity of the tube.
The conduction heat transfer coefficient presented in water/ice U cd,w [W/m 2 °C] is calculated as The overall heat transfer coefficient U [W/m 2 °C] is obtained by considering the convection and conduction heat transfer coefficients ( 17)-( 19) and the electric analogy presented in [38] (as shown in Figure 6), yielding where β is a correction factor.In this paper, β = 1.5, and this value was obtained through a heuristic process as in [30] to fit simulation with experimental results.
It is important to emphasise that the energy balance and the equations for heat transfer coefficients presented in this section consider that the thermophysical properties of the HTF and of water/ice are temperature dependent.

| Model discretisation
System components such as heat exchangers or TES tanks, where heat transfer occurs due to a HTF flowing through a channel, have been previously modelled in the literature using a stratification or discretisation approach [39][40][41][42].For the icebased tank, a division into nodes captures how the heat varies through the tubes during the heat transfer process due to the temperature changes of the HTF and of water/ice.For a pair of tubes, this means that the temperature of the HTF and of the water/ice closer to the HTF inlet are modified faster than nearer to the outlet.These dynamic variations are well described by dividing the length of the tubes into nodes.In the discretisation method, a number of control volumes are specified to, in turn, define the boundaries of the heat transfer process as described by ( 6), ( 9), (11), and (12).A schematic description of this method when applied to a single spiral tube in the ice tank is shown in Figure 7.
As discussed in Section 2, having two tubes per horizontal level requires conducting energy balance for each tube.Since the tank considers several horizontal layers, the number of equations required for energy balance is significant.However, as the flow for a pair of tubes is assumed the same as in any other pair of tubes throughout the tank's height, the heat transfer process occurs under equal conditions for all pairs of tubes.To explain this, Figure 8 shows the two unrolled tubes for a single horizontal level during a charging process.The cold fluid input in the opposite extremes of the tubes freezes the water around them starting from the closest node, as shown in Figure 8(a).The node numbering is shown in Figure 8(b).Tubes a and b are split into the same number of nodes and the heat transfer process occurs in similar conditions but in opposite directions.Hence, for simplicity, the energy balance considers only the control volumes of the HTF and water/ice for a single horizontal level comprised of tubes a and b.
By applying the discretisation method, a set of nonlinear differential equations for energy balance is obtained, with four equations per node and N denoting the total number of nodes.The 4N-th order system can be described in matrix form to adopt a state-space notation, where temperature is the state variable for each equation.Thus, the mathematical model is described by (21), where the input temperature to each control volume is the HTF temperature from the previous node (T f,k−1 ) for tube a and the next node (T f,k+1 ) for tube b.This means that the contribution of the heat source (i.e. the HTF) is affected by the thermal conditions of the preceding node in the direction of the HTF.
The energy storage capacity of an ice-based TES tank is given by the amount of water/ice and its LHV.The total energy E tot stored when the tank is completely charged is defined by where m w [kg] is the total mass of water and ΔH L,m [ J/kg] is the LHV of water/ice (for melting-solidification).The LHV is the total heat released or absorbed during the phase change of a material.This internal thermal energy is described by the enthalpy curve.Two phase changes of interest are presented in water.These are shown in Figure 9, where the LHV for meltingsolidification is denoted by ΔH L,m .Conversely, the LHV for evaporation-condensation is denoted by ΔH L,e .A meltingsolidification LHV suitable for the dynamic model of the ice tank presented in Section 3 is required, which will, in turn, enable the calculation of the internal energy of the water/ice.This is because the energy balance described by ( 9) and ( 12) is specified in terms of specific heat as opposed to enthalpy.The following relationship between specific heat and enthalpy exists: where H w is the enthalpy of water, and T 1 and T 2 are the temperature limits where H w is calculated (given a curve for specific heat of water c p,w ).From (23), it is possible to determine a specific heat curve by differentiating an enthalpy curve with respect to the change of temperature.This curve defines the LHV of water/ice F I G U R E 6 Top and frontal cross sections of the control volumes of tubes a and b.Blue and red arrows show the radial direction of the heat.For a charging process, blue arrows denote heat flow towards the HTF, whereas for a discharging process, red arrows denote heat flow towards the water/ice.The white and black arrows show the heat transfer between the ice/water control volumes of the tubes.An electrical analogy is shown at the bottom, where U w , U cv , U cd,t , and U cd,w are the heat transfer coefficients defined in ( 8) and ( 17)-( 19) during phase transition in the neighbourhood of the melting temperature (i.e.0°C).However, accurate curves are difficult to determine as there are instances in time when solid and liquid phases coexist.A liquid fraction analysis would determine the liquid and solid portions during a melting or solidification process so as to define an accurate enthalpy curve [7]; however, such an exercise falls out of the scope of this paper.Instead, specific heat curves for discharging and charging processes are defined by assuming that the temperature of the control volume of water/ice directly determines the value of specific heat without considering the liquid fraction [5,30,43].The derived curves are modelled using a probability density function (PDF), which allows to shape the specific heat curve by adjusting a few parameters.This approach has been implemented with good results in [30,36].The temperature is the independent variable in the PDF, and this is used to calculate the LHV.The specific heat curves are defined in terms of the PDF φ T ð Þ as Equation ( 24) has been adapted from the general equation of the specific heat proposed in [30] to consider constant liquid and solid states.Coefficients a 0 , a 1 , a 2 , σ, τ and γ in (24) and ( 25) are dimensionless parameters which have been heuristically adjusted to replicate experimental results with the presented mathematical model.The values used in this paper are shown in Table 1.
Figure 10 shows the specific heat curves derived in this paper, together with their corresponding enthalpy curves.The specific heat curves enable to calculate the amount of thermal energy that has been released or absorbed during the phase transition by knowing the temperature of water/ice.

| SIMULATION RESULTS AND DISCUSSION
The mathematical model presented in Section 3 was verified using experimental results available in [33,34], which have been obtained with the ice-based tank model ICEBANK 1098C [35].Parameters for a single tube are provided in Table 2.A 34% mass water-glycol mixture is used as the HTF (see Appendix A for further details).
A special tank setup is used in [33,34] to carry out the experiments, where the top 16 levels out of the total 34 are blocked and the HTF circulates through the bottom 18 levels only, which are submerged in water.The total volume of water of the tank is V w = 3710 L [35] and the partial mass of water/ ice is approximately m w = 1876 kg; this implies an effective volume of V w,e = 1876 L = 1.876 m 3 (nearly half of the volumetric capacity) [33].The storage capacity for such a configuration is 174.05 kWh, which is almost half the 350 kWh established by the manufacturer.To avoid discrepancies with respect to the experimental results, the parameters provided in Table 3 are adopted in this paper.These have been determined using the procedure included in Appendix B.
T A B L E 1 Parameters used in (24) and (25)   The ice-based tank model has been implemented in MATLAB/Simulink using S-functions and look-up tables.These encapsulate the temperature dependence of the thermophysical properties of water/ice and the HTF.As discussed previously, the discretised model consists of 4N ODEs.For this purpose, the system described by ( 21) is coded employing nested functions due to the link between the temperature of the HTF at any node with the temperature of the HTF at a previous or subsequent node.Using a state-space notation (where temperature T is replaced with x), MATLAB's ODE solver is employed.
The pseudo-code in (26) represents the mathematical model of the tank as a set of nonlinear ODEs, with solutions stored in array x.Such an approach enables the implementation of a system with an arbitrary number of nodes, given by N. The first and last nodes are defined separately as their input temperatures are the overall input temperature of the system T f,in (and not the output temperature from a previous node for tube a or from a subsequent node for tube b).Odd entries of array x provide node temperatures of the HTF, whereas even entries the temperatures of water/ice.Simulation results for a discharging process are compared with the experimental results reported in [33] in Figure 11.The number of nodes has been selected as N = 20, which implies that the ice-based tank is defined by 80 ODEs.In the figure, subscript 'e' is used to denote experimental results, 'in' for input, 'o' for output, '20,a' for the last node of tube a and '1,b' for the first node of tube b in the simulation results, 'f ' for HTF, and 'w' for water/ice.The conditions of the first experiment are a mass flow rate and input temperature of the HTF of _ m f ¼ 3689 kg/h = 1.027 kg/s and T f,in,e = 12.8°C.
As it can be observed in Figure 11, a good agreement between experimental data (T f,o,e ) and simulation results (T f,20,a and T f,1,b ) is achieved.Since a similar mass flow rate of HTF occurs in both tubes a and b, the mixing temperature of the flows is their average temperature.However, as the output temperatures are equal, the mixing temperature is thus not shown.It is important to highlight that the small notch exhibited by the input temperature T f,in,e in the experimental results from [33] is reproduced by the simulation results.
The average cooling power provided in [33] is also compared.This is calculated with Screenshots showing the model implementation in MATLAB/Simulink are provided in Figure 12.The input  temperature of the HTF and the experimental data (output temperature and power) are included in the simulation model using look-up tables.Equation ( 27) is solved during the simulation.
A flowchart describing the use of equations ( 14)-( 21), ( 24), ( 25) and ( 27) is presented in Figure 13.All equations are computed every time step.Thus, the vector ODE ( 21) is solved and the solution is used to calculate the thermal power using (27).
The graph shown in Figure 14 shows a comparison between simulation and experimental results of the cooling power for the same experiment as in Figure 11.The cooling power as reported in [33] (P e,a ) presents a mismatch with respect to the one calculated from simulation results (denoted by P s ).Nevertheless, if cooling power is calculated directly with equation ( 27) using the experimental values of T f,in,e and T f,o,e instead of T f,o,s (yielding the trace denoted P e,b ), a better agreement with the simulation results is achieved.As it can be concluded, the model presented in Sections 3 and 4 is capable of reproducing the experimental behaviour.
The sequence of the ice melting process can be seen in detail through the temperature of each node.Figure 15 shows the HTF and water/ice temperatures of all nodes for tube a during the discharging process.The node temperature varies gradually depending on the specific location of the node.As expected, the HTF and the water/ice of the first nodes increase their temperature faster than those nodes located at the end of the tube.This difference is given by the higher temperature of the HTF at the beginning of the tube.The same thermal performance is exhibited in the nodes of tube b but with an opposite pattern (not shown).
Note: The rationale behind the choice of 20 nodes is that incrementing this number would not produce a significant improvement in the accuracy of the simulation results.
Figure 16 shows the simulation results of the discharging process using different number of nodes (5,10,20,30,50,100).For clarity, only results for tube a are shown.Subscript 'xN' in the figure stands for the number of nodes considered in the model.As it can be observed, the difference between the output temperature of the HTF of the model with 20 nodes (red traces) with respect to a model with 100 nodes (dashed grey traces) is about 0.1°C.Comparing the water/ice temperature in the last node of tube a for models with 20 and 100 nodes shows a marginal difference of 0.03°C.The value of 20 nodes is thus adopted as a compromise between model complexity (i.e.80 ODEs) and accuracy.
Although a unique discharging process recording temperature data is presented in [33], the cooling power graphs for four additional discharging runs without any temperature results are included in the reference.Table 4 shows the operating conditions for these additional experiments.Simulations were conducted to further verify the mathematical model presented in Sections 3 and 4 using the experimental results reported in [33].For all simulations, identical initial conditions of −5°C were adopted for all node temperatures (both HTF and water/ice).
Figure 17 shows the results of Experiment 1 according to the conditions in Table 4.The graph in Figure 17(a) shows the input temperature for the experiments T f,in,e and the simulated temperatures of the HTF and water/ice at the final node of tube a (T f,20,a and T w,20,a ) and at the initial node of tube b (T f,1, b and T w,1,b ).Similar responses of the HTF and water/ice temperatures to those of the previous simulation were obtained.The cooling power is shown in Figure 17(b), with subscript '1' denoting Experiment 1.It can be seen that the power curve obtained using the simulation results and equation (27) agrees on well the experimental data (P e,1 ).
Figure 18 shows a comparison between the cooling power P e obtained in Experiments 2-4 as reported in [33] and the cooling power obtained with the simulation model presented in this paper.As in Experiment 1, temperature information has 4 Comparison between experimental and simulated mean cooling power.The small notch presented in P e,a at 110 min comes from the experimental results in [33]  For completeness, the charging process as reported in [34], incidentally adopting the same tank setup as in the discharging experiments, was replicated here with the simulation model.A vapour-compression chiller system supplies the cold HTF.The operation conditions are T f,in,e = − 3.7°C and _ m f ¼ 0:6817 kg/s and the initial condition of temperature of the HTF and water/ice in all nodes is 10.5°C.The inlet temperature profile is shown in Figure 19(a), where it can be noticed that the temperature provided is not constant throughout the duration of the charging process.It should be highlighted that the time required to charge the ice-based tank in [34] is large (nearly 25 h).These operating conditions do not necessarily represent the real operation of a cooling centre where the tank would be charged at night.However, the comparison of temperature profiles shows that the simulated HTF temperature of the last node of tube a (T f, 20  The experimental input temperature (T f,in,e ) and the temperature of nodes 1, 10 and 20 of tube a are shown in Figure 19(b), where subscript 'a' has been dropped from the labels for clarity.As it can be seen, the HTF temperature in the first node T f,1 decreases faster than in any other node-due to the lower temperature of the HTF at the beginning of the tube.As the temperature of the HTF increases through the tube's length, the temperature in the last nodes decreases more slowly than in the first nodes.A similar thermal performance occurs in the nodes of tube b but with an opposite pattern (not shown).The large decrement in the input and output temperatures presented after 3.3 h could be attributed to a subcooling phenomenon or to the temperature delivered by the cooling source (vapour-compression chiller).
The subcooling phenomenon that may be exhibited during charging processes provokes a sudden decrement of water temperature when the ice starts to build [7].This directly affects the temperature of the HTF, and is observed in Figure 19 at 5 h.The effects of subcooling are challenging to represent with a 1-D model and, together with broken ice movement towards the top of the tank during discharging processes, fall out of the scope of this paper.
Note: The mathematical model of the ice-based tank presented in this paper may be implemented in any software engine with an ODE solver.The model defined by (21) can be built based on the pseudo-code given by ( 26) and equation (27), the dynamic calculation of the overall heat transfer coefficient, the specific heat curve of water/ice defined by ( 24) and ( 25), the temperature-dependent thermophysical properties of the HTF provided in Appendix A, and the input parameters defined by Tables 2 and 3.The model should update the input temperature and the mass flow rate of the HTF during the simulation time to consider any changes in operating conditions.

| On the computation time and time steps
The solution of a set of ODEs in an engine solver requires establishing a time step.This is a critical parameter as the  4 accuracy and speed of the solution will be influenced by this value.As a result, the computation time will be also affected.
The simulation results presented in Section 5 were obtained using a variable time step-which is available in the MATLAB/Simulink solver.The rationale behind this choice is that small time steps are not required as the thermal dynamics of the ice-based tank are slow.However, for completeness, the effect of this simulation parameter over the presented model was also investigated.
The discretised model with 20 nodes was simulated with a variable time step and with fixed time steps of 0.1 and 1 s.The output temperature of the HTF of tube a is shown in Figure 20, where subscripts '0.1' and '1' stand for fixed time steps of 0.1 and 1 s and 'v' for a variable step.Although the accuracy of the simulation results is not affected by the time step selection (see the zoomed-in area between 135 and 155 min), the computation time (i.e.run time) varies dramatically.This information is summarised in Table 5.
To further examine the computation time of the presented model, simulations of a discharging process considering different numbers of discretisation nodes have been conducted.These simulations considered the same conditions as for the simulations reported in Figure 16, which were carried out using a variable time step.In this case, a fixed time step of 1 s was adopted.This value of time step was chosen for consistency with the information in Table 5.The run times for all simulations are summarised in Table 6.The maximum and minimum computation times with a variable step were 81.09 s for 100 nodes and 4.78 s for 5 nodes.All simulations with a fixed time step of 1 s had a computation time of around 90 s irrespective of the number of nodes used.
Comparing the simulation run times given in Table 6, corresponding to an 8.33-h discharging process of a whole TES tank, with the 5 min required to simulate a nearly 2-h charging process of just a single small PCM sphere with a 3-D model [45], the computation times afforded by the 1-D model here presented are deemed appropriate.
The information summarised in Tables 5 and 6 justifies the selection of a variable step to conduct the simulations in Section 5.The central processing unit (CPU) used for all simulations was an Intel Core i7-10,610U CPU @ 180 GHz.

| Error analysis
For completeness, an error analysis was carried out between experimental results reported in [33,34] and the simulation results obtained with the presented model.Given that the output temperatures of the last node of tube a and the first node of tube b are the same, the analysis is restricted to tube a only.
Although maximum error values of 2.66°C and 10.87 kW (see Table 7) might indicate large discrepancies between simulation results and experimental data, these errors are exhibited only at the beginning of the processes or when subcooling may be taking place (which is a phenomenon not accounted by the presented model).The small values of MSE F I G U R E 2 0 Simulations of the discharging process (output temperature of the HTF, last node of tube a) using three different time steps: fixed step of 1 s (T f,1 ), fixed step of 0.1 s (T f,0.1 ) and variable step (T f,v ) T A B L E 5 Effect of different time steps in the simulation run time for a discharging process

Simulation step
Run time (sec) Step size (sec) demonstrate good agreement between simulation and experimental results.

Note:
The initial values of power for the plots in [33] were not considered in the error analysis as they were set to 0 kW in the reference.This may be unrealistic since a maximum difference between the input and output temperature of the HTF would exist by the start of the tests.To avoid a misleading interpretation of the errors, these have been calculated from the moment when power has a maximum value in the plots shown in [33].

| ON THE LIMITATIONS OF THE PRESENTED MODEL
Although an accurate description of the thermal behaviour of an ice tank is achieved with the model presented in this paper, it is important to acknowledge its limitations.Firstly, the fractions of liquid and solid water are not described by the model.Instead, a specific heat curve is used to represent the phase change transition.
Momentum and mass balances are not considered when mathematically formulating the model.Therefore, mechanical conditions such as pressure increment during ice growth or density changes in the HTF are not calculated.
Due to the adoption of lumped control volumes for energy balance, the modelling approach has a limited heat transfer calculation as opposed to models that do not neglect spatial components in the y and z directions.Therefore, a correction factor of the overall heat transfer coefficient is required to improve the heat transfer analysis.
Verification of the model with a complete experimental charging setup of the tank would provide further confidence in the modelling approach adopted in the paper.This can be complemented with an estimation of the state-of-charge of the tank as in [46].

| CONCLUSIONS
In this paper, a 1-D dynamic model of an ice-based TES tank has been developed.The mathematical model has been verified against experimental data from a commercial configuration available in the literature.The model is sufficiently detailed to accurately describe the heat transfer phenomena in TES systems, but significantly simpler than 2-D models and the computationally demanding 3-D models.The simplicity of the developed model is based on the assumption of an equal temperature distribution and equal hydraulic conditions for all pair of tubes inside the tank, which reduces the analysis of the heat transfer phenomena to a single pair of internal tubes only instead of 68 tubes.
The model captures the temperature dependence of the thermophysical properties of the HTF and of water/iceespecially the specific heat and enthalpy.The PDFs used to dynamically describe these properties offer a powerful solution to simulate the phase change.Since the specific heat or enthalpy curves of commercial PCMs are not commonly available, the presented modelling approach relieves this key shortcoming.This, in turn, enables the model to be adapted to different types of PCM.
The 1-D ice-based TES tank model has been built in MATLAB/Simulink.Its implementation is simpler compared to 3-D models as it is described by ODEs instead of partial

Maximum error MSE
Output temperature of HTF, e T,a differential equations.An advantage of the model presented in the paper is its portability and versatility: the ease of its implementation as a set of nested ODEs may be easily replicated in any simulation software incorporating an ODE solver.Although the model is compact, it is sufficiently comprehensive to be adapted to different internal geometries of a TES system and of different PCMs.By suitably modifying parameters in the mathematical formulation, the model can be extended from cooling to heating configurations.The portability of the model is suitable for system level studies of cooling energy centres or advanced industrial processes.

APPENDIX A
The thermophysical properties water-glycol at 34% used to solve (21) are obtained from [47].These are provided in polynomial form as follows: where T is the temperature of the fluid in °C.

APPENDIX B
The volume of water/ice used in this work is determined from the dimensions of the control volume defined in Figure 4 and the length of the tubes.According to [33,34], half of the volume of the tank is about 1.876 m 3 .However, the experiments were carried out with 18 levels out of the total 34.The proportional volume of water with respect to the total volume (V w = 3.710 m 3 ) is thus 1.9641 m 3 .The external and internal radii of the control volume are estimated considering the approximated tube dimensions in [33].The external radius (r ice = 23.939 5 mm) is defined considering a distance of 60 mm between the centres of the adjacent rolled tubes.The total value of the tubing length in each horizontal level (considering the two adjacent tubes per level) is established as 65.1163 m, which is close to the approximate value of 70 m defined in [33].
Given that there two tubes for each horizontal level, the length of each tube is L = 32.558 1 m.This way, the total volume of water/ice for each tube is given by where r o is the external radius of the tube.The total volume of the control volume of water/ice is obtained with V T,ice = 2V ice N L = 1.878 3 m 3 , where N L = 18 is the number of tube levels, and two single tubes are considered for each level.This volume is close to the proportional water volume for 18 levels.A slight increase in the volume of water/ice leads to a small increment in the storage capacity of the tank (E tot = 174.264kWh).This was computed with ρ w = 1000 kg/ m 3 .This value is very close to the amount of energy of 174.05 kWh defined for the ice-based tank setup in [33]. 18-BASTIDA ET AL.

1
Schematic of a cooling energy centre with two ice-based TES tanks U R E 2 (a) Internal geometry of a single ice-based storage tank showing the direction of the fluid in each level.(b) Location of the two external and two internal headers [31] F I G U R E 3 Top view of the ice tank during discharging (left) and charging (right) processes.The white shading denotes ice and the cyan shading represents water BASTIDA ET AL.
where _ E f ;b is the rate of change of energy of the HTF for tube b, _ E w;b is the rate of change of the thermal energy of water/ice for tube b, and the heat loss _ E l;b for tube b is defined by

9 F
Enthalpy of water showing two phase changes and LHVs F I G U R E 7 Control volume (nodes) boundaries of one tube I G U R E 8 (a) Unrolled tubes with fluid directions and ice built in the first nodes for tube a and the last nodes for tube b, where the white shading represents ice formation.(b) Control volumes through the length of the tubes

F
I G U R E 1 0 (a) Specific heat curves derived in this paper.(b) Enthalpy curves obtained from the integration of the specific heat curves.Subscripts 'dch' and 'ch' stand for discharging and charging processes BASTIDA ET AL.

F I G U R E 1 1
Comparison between T f,o,e and the temperatures of the last node of the HTF of tube a (T f,20,a ) and of the first node in tube b (T f,1,b ).T f,in,e , T w,20,a and T w,1,b are also plotted.The small notch presented in T f,in,e at 110 min comes from the experimental results in [33] F I G U R E 1 2 Screenshot of the model implementation in MATLAB/ Simulink: system inputs (top), S-function of the 1-D dynamic model of the TES system (middle), and calculation of power (bottom) BASTIDA ET AL.

5
Simulation results of tube a. Node temperatures during discharging process: (a) fluid temperatures; (b) water/ice temperatures.The small notch in T f,in,e is reproduced by the temperature of HTF nodes F I G U R E 1 3 Flowchart describing the mathematical relations of the 1-D model of the TES tank been omitted.Although there is a slight mismatch between simulation and experimental results, this is deemed acceptable considering the lack of data of input and output temperatures of the HTF.
,a ) and the first node of tube b (T f,1,b ) with the output temperature of the experiment (T f,o,e ) agree on well.

6 T A B L E 4 7
Comparison between ice-based TES tank models with different number of nodes.(a) Output temperature of the HTF.(b) Water/ ice temperature in the last node of tube a Operating conditions of discharging processes Discharging process for T f,in,e = 14.6°C and _ m f ¼ 1:131 kg/s.(a) Input temperature of the HTF and output temperatures of the HTF and water/ice.(b) Comparison between experimental and simulated average cooling power BASTIDA ET AL.
U R E 1 9 (a) Experimental and simulated results of the charging process.(b) Temperature of fluid and water/ice of 1, 10 and 20 nodes of tube a during a charging process F I G U R E 1 8 Comparison of the cooling power for the operating conditions in Table Figure 21(a) and Figure 21(b).The errors between the power curves in Figure 17(b) and Figure 18 (e s,x = P e,x − P s,x ) are shown in Figure 21(c).Finally, the error of the output temperature of the HTF for a charging process (e T,ch = T f,o,e − T f,20,a ) is shown in Figure 21(d).The maximum values of error and the mean squared errors (MSE) are summarised in Table7.Although maximum error values of 2.66°C and 10.87 kW (see Table7) might indicate large discrepancies between simulation results and experimental data, these errors are exhibited only at the beginning of the processes or when subcooling may be taking place (which is a phenomenon not accounted by the presented model).The small values of MSE

1
Error between experimental results (from[33,34]) and simulation results with the presented model: (a) output temperature of the HTF (tube a) for a discharging process, (b)-(c) calculated power during discharging processes; and (d) output temperature of the HTF (tube a) for a charging processT A B L E 7 Maximum value of errors and MSEs between experimental and simulation results

T A B L E 2
Parameters of ice tank (single tube) Note: The model does not calculate the liquid or solid fractions of water during phase change.Thus, an average value of k w is assumed during the phase transition, where 2.1 W/m°C is used for the solid phase and 0.6 W/m°C for the liquid phase.