Numerical Modeling of Thermal Storage Performance of Encapsulated PCM Particles in an Unstructured Packed Bed

: In this work, a numerical model to study phase transition in unstructured packed beds of phase change particles for latent heat thermal energy storage is presented and implemented into the open-source CFDEM R (cid:13) coupling program. The governing equations for both ﬂuid and solid particles are discussed. The presented model is validated using experimental data. The validated model is used to study the effect of bed structure (e.g., using multiple sizes of particles and mixing of multiple conﬁgurations of particles) on the charging and discharging of PCM bed. It is found that using smaller particles leads to faster charging and discharging of the bed. A decrease of 26% was achieved by changing half the bed into particles with half the diameter. Furthermore, it is observed that placing small particles downstream has more effect on increased charging speed than placing them upstream when charging the bed fully. A decrease of 20% in charging time and a decrease of 13% in discharging time were observed in the conﬁgurations that were tested.


Introduction
The share of renewable energy in total primary energy supply is on the rise with the need to move away from the dependency on fossil fuels and to reduce CO 2 emissions. To achieve the target of the Paris Agreement, CO 2 emissions need to be reduced significantly [1]. The total energy usage of the residential and service sector accounts for 35.5% of the total global energy consumption, from which 75% is used for space and domestic heating [2]. Solar energy can provide heating without the need for fossil fuels by using solar boilers. However, the major problem of solar energy is the gap between peak power production and the peak power demand.
To bridge this gap, there is an urgent need for introduction of heat storage systems. Two common forms of heat storage are Sensible Heat Storage (SHS) and Latent Heat Storage (LHS). SHS uses sensible heat to store thermal energy, meaning that the energy stored in the Thermal Energy Storage (TES) is simply used to heat up a material from which this heat can later be extracted. This version of TES is cheaper and simpler than LHS, which uses both sensible heat and latent heat [3]. When LHS is used, the storage medium undergoes phase change in the temperature range in which the storage system operates. The advantage of this is that, for similar storage capacity, LHS offers a higher thermal energy storage density compared to SHS. Another advantage is that during phase change the temperature is constant, thereby allowing a constant temperature difference between heat transfer fluid (HTF) and energy storage material.
Various LHS systems exist which employ phase change material (PCM) in a variety of ways [4][5][6]. Of these, a packed bed of PCM particles provides the highest heat transfer area and thus allows for both quick charging and discharging. These LHS systems are of interest for use in both large-scale concentrated solar power (CSP) plants and small-scale domestic applications. In the case of CSP, the TES is a way to store high temperature heat used to drive a generator to produce electricity. Mohammadnejad et al. [7] studied the use of different flow velocities and bed porosity during discharging and the effect of different configurations of particles. Senthil et al. [8] studied the positioning of the rectangular solar rectangular solar receiver in a CSP setup.
In domestic applications, on the other hand, TES is used to store heat purely for hot water consumption in the form of domestic heating and sanitary applications. For example, Pakrouh et al. [9] studied a system that uses water as HTF and paraffin wax as PCM. This form of domestic energy storage is often done by simple SHS systems. However, the use of PCM can increase the efficiency of the TES by storing the energy not just in sensible heat but also in the latent heat of the PCM.
To help in the design of the LHS, multiple experimental studies were conducted, for instance, by Barrientos et al. and Nullasamy et al. [10,11]. Conducting such experiments is often costly, while the use of numerical models can offer a cheaper method to provide design criteria for LHS system, limiting the number of experiments needed. Another advantage is the amount of detailed information that can be extracted from the model. Furthermore, modeling allows direct control over all variables within a system without a concern for deviations and measurement accuracy. Besides, modeling an unstructured packed bed gives the possibility to measure the temperature of every single particles, which is practically unfeasible in experimental measurements.
Various methods can be used to model a beds of phase changing particles (e.g., single-phase modeling, two-phase modeling, and continuous or discreet solid phase), as shown in the overview by Garcia et al. [12]. Mao and Zhang [13] combined three PCMs for use as a heat storage medium. They studied numerically the impact of of particle diameter, porosity, and height-to-diameter ratio of the storage tank on the total storage energy, storage capacity ratio, axial temperature curve, and utilization ratio of the PCM. Khor et al. [14] investigated both numerically and experimentally the effect of the fillers in the void space between the macro-encapsulated PCM of a packed bed to increase the storage capacity of the tank. The PCM particles are placed in the tank in a structured way. They found void space reduction with the granular materials allows maintaining the storage efficiency while achieving an increase in phase change material encapsulation size to reduce its overall cost. Duan et al. [15] used a numerical model to find the suitable structured packing for a packed bed of PCM capsule. They found that FCC arrangement (face center cubic packing) has superior characteristics, such as high charging rate, low investment cost, large total heat storage energy, and stable output temperature, compared with BCC arrangement (body center cubic packing) or simple cubic packing.
The response time in a heat storage system is an essential parameters that should be as high as possible. Therefore, the system must be designed in such a way to enhance the heat transfer between heat transfer fluid and PCM particles. The bed structure has a significant impact on the heat transfer through a packed bed. Therefore, the effect of this parameter must be well studied to find the optimum structure for a faster response time and a higher energy density. However, there is a gap in the literature to address these challenges. The objective of this work is to investigate the effect of the bed structure on the charging and discharging of a packed bed of PCM. For this purpose CFDEM R coupling, from the open-source CFDEM R project [16], is used to calculate the heat transfer between particles (i.e., PCM) and the surrounding fluid phase. Moreover, the model has been further developed to calculate the phase transition within the particle. The outcome of this study provides detailed information that can be used by the designer to further optimize the PCM storage tank.

The Mathematical Model
To model an unstructured bed of PCM particles, the model first proposed by Schumann is used [17]. This model is a two-phase model meaning that the particle and fluid phases are treated by separate equations which are coupled. Within the present model, a packed bed is considered as an ensemble of a finite number of particles (i.e., discreet phase). Transient conservation equations for energy is solved for PCM particles individually. Applying this model to all particles of a packed bed forms the entire packed bed process as a sum of the individual particle processes. The flow through the void space in the bed is modeled as a flow passing through a porous media (i.e., CFD domain), while the interaction between the solid and the fluid phases by heat and momentum transfer is taken to account. To avoid porosity equal to zero in the CFD domain, particles must be smaller that the CFD cells. Water is used as the HTF in all cases and is considered as an incompressible fluid and continuous phase; the governing equations therefore can be written as follows.
Continuity equation for the fluid: Momentum equation for the fluid: where is the void fraction, t is the time, u f is the fluid velocity vector, p is the pressure, ρ f is the fluid density, and τ is the stress tensor. This particular form of the momentum equation is based on the so-called Model III as given by Zhou et al. [18]. The last term on the right hand side is the coupling term of momentum transfer between the particles and the fluid with: where u p is the particle velocity vector and |∑ i F d | is the sum of all the drag forces in the computational cell with volume V cell . Energy equation for the fluid: where T f is the fluid temperature, c p, f is the fluid heat specific heat capacity, and k f is the thermal conduction coefficient of the fluid. The last term on the right-hand side is the heat transfer between the particles and the fluid due to convection. Radiation from the particles are neglected, thus it is assumed that convection is the only mode of heat transfer between particles and fluid. The particles have their separate equations to govern conservation of momentum and energy.
Momentum equations for the particles: where m i is the mass of the ith particle;ẍ i is the acceleration vector of the particle; F i,n and F i,t are, respectively, the force in normal direction and the force in tangential direction due to particle-particle interaction; F i,d is the drag force on the particle; g is the gravitational acceleration vector; ρ p is the particle density; and V p is the particle volume.
where I i is the moment of inertia of the particle, ω i is the angular momentum vector, and r i is the moment arm vector. Energy equation particles: where c p,p is the particle heat capacity, T p is the particle temperature, and q p f is the heat source/sink therm. It is assumed that the particles act as a lumped system and therefore have a single temperature. Furthermore, the particles are frozen in place during fluid-particle simulations and therefore the bed is a fixed bed. Thus, all forces in Equations (5) and (6) are set to zero during particle-fluid simulations. However, the drag force is still calculated and taken into account in the fluid momentum equation. Freezing the particles in place is a valid assumption since in reality their movement would be limited.
It is also assumed that the only mode of heat transfer which the particles experience is convection with the fluid. Particle-to-particle conduction and radiation are neglected, as these are small compared to the convection term. Lastly, it is assumed that the walls of the container are well insulated so losses to the environment are neglected as well.
The heat transfer coefficient between particles and fluid can be determined using Nusselt correlation from Li and Mason [19] with coefficient variable n being 3.5: where the Reynolds number is and the Prandtl number is Through Equation (8), the heat transfer coefficient U p f can be determined which can then be used to calculate the specific heat flux between a particle and the fluid due to convection.
To take into account phase change, the enthalpy method is used by taking the heat capacity of the particles as a function of temperature [20,21]. To ease the modeling of this transition, the heat capacity curve has to be approximated; for this, multiple models exist [22], as shown in Figure 1.
In the validation with measurements [11], the exact profile of the Cp as a function of the temperature for the used PCM is not known. However, it was found that the right triangle method was the most agreeable with experimental results. Since the same PCM is used for validation as is used in the models presented here, the right triangle method is used, leading to the heat capacity as follows: This change in c p has been incorporated into the CFDEM R coupling program, thus allowing for the modeling of phase transition in a packed bed of particles.

Validation
The model was validated using the data available from experiments carried out by Nallusamy et al., which concerns a structured simple cubed bed filled with Paraffin wax particles [11]. The results of the model compared to the experiment are given in Figures 2 and 3.
In these figures, it can be seen that their general behavior is the same, and they are in a relatively good agreement. The particle heating curve in Figure 2 shows that initially the particle temperature is overestimated by the model. This is due to the experiment measuring with temperature sensors inside the particle. Since the model assumes a lumped system, the temperatures of modeled particles are overestimated compared to the interior of the actual particle. The particle centers experience a delay in the heating up due to the thermal resistance.
After a while, the numerical curve flattens and the phase change start in the model. Since the heating curve in the experiment is also not straight at this point, we can conclude that in actuality there is also some phase change taking place. However, the model seems to overestimate the amount of latent heat storage occurring between about 320 and 330 K, which is due to the c p (T) curve being only an approximation. This means that effect of phase change at lower temperatures is exaggerated in the model. This can also be seen at 334 K, where the experimental temperature line is almost flat while the numerical line is climbing ; this suggests that a greater part of the phase change occurs at higher temperatures in the experimental case.
This effect of exaggerated phase change at lower temperatures can also be seen in the temperature of the HTF in Figure 3. The increase of latent heat storage at lower temperatures causes a greater amount of heat transfer between fluid and particle. This is because a greater difference in temperature between particles and HTF is maintained for a longer time. Therefore, the numerical model temperature is higher. The flattening effect on the temperature curve of the HTF due to phase change is also spread out more due to the overestimation of latent heat storage at lower temperatures.
From this, one can conclude that the model is still quite accurate; however, the effect of the approximation of the c p (T) curve and the lumped system assumption can be seen. Using a more accurate c p (T) curve, which requires dedicated measurements, could alleviate this difference. Nevertheless, if the effects of the approximation of specific heat capacity c p (T) and the lumped system assumption are kept in the back of ones mind, the general behavior of the system can still be modeled accurately and conclusions about the systems behavior can be drawn.

Model Setup
The particles are made of paraffin wax with the thermodynamic properties as given in Table 1. All particles are assumed to be spherical with a diameter of 55 (for large particles) or 27.5 mm (for small particles). Water is used as the HTF with the inlet temperature of 343.15 K during the charging phase and 293.15 K during discharging. The initial temperatures of the particles are 293.15 and 343.15 K for charging and discharging cases, respectively.
In all cases, the flow rate is constant at 2 L per minute at the inlet and move from bottom to the top. The energy storage required is based on heat demand of a single household, which is approximated to be 120 L of water at 60 • C or 333.15 K. By considering an overcapacity of 20% and the stored energy in the PCM as well as the sensible heat stored in the water around the particles, 670 large (55 mm) particles are required to meet storage demand.
An overview of the different configurations is given in Figure 4. The unstructured beds are generated using DEM (Discrete Element Method) by free falling of particles. In all cases, the total volume and consequently the total mass of the PCM are the same. During all simulations, it is assumed that the wall is perfectly insulated.  For the CFD domain, a tank of 0.5-m diameter and 11-m height is defined. To store all the particles, a smaller tank would have been sufficient but adding extra CFD cells before and after the bed increases the stability of the model. The combined CFD and DEM domain is shown in Figure 5. Here, one can see half of the CFD domain with the DEM particles. The flow goes from the inlet at the bottom to the outlet at the top for charging (heating up) cases and the reverse for discharging (cooling) cases.

Effect of Particle Size on Performance
To find the effect of particle size on TES, Configurations I-IV are compared for a charging cycle. During the charging process, HTF is pumped from the bottom through the packed bed. The average temperature of particles during the charging is giving in Figure 6. The average temperature of the particles is weighted by both the mass and c p,p value of the particles, thus Figure 6 shows that the average temperature of the particles increases faster when there are more small particles. At 3.5 h, Configurations III and IV have a temperature of about 340 K while the particles in Configuration II have a temperature of 338 K and Configuration I has a temperature of 337 K. This is due to the increase in the area available for heat transfer which provides faster heat transport; Configurations III and IV have similar temperatures because they have the same surface area. Here, it can be seen clearly that the small particles heat up significantly more quickly. This also explains why in Configuration I the average temperature shows a more clear flattening compared to Configurations II and III. In fact, the addition of small particles results in faster heating up. This compensates for the large particles which are still in phase transition and therefore their temperature does not increase. In Configuration IVa, this effect is less present since the large and small particles are not intermixed. The energy stored in the PCM, which is the sum of all sensible and latent heat inside the particles, can be seen in Figure 7. In this graph, it can be seen that the energy contained within the particles also increases more quickly as a result of more heat transfer area due to the addition of the small particles. This indicates that, to increase the charging speed, the particle size should be as small as possible. To achieve a fully charged bed, Configuration I requires a charge time of 6 h 17 min, Configuration II requires 6 h 5 min, Configuration III requires 5 h 17 min, and Configuration IV requires 4 h 38 min. Thus, for the same mass of PCM, Configuration IV charges more than 1.5 h more quickly than the others. It should be noted that Configurations III and IV have the same surface area available for the heat transfer. The difference between the energy stored within the bed is caused by the different placement of the smaller particles. Therefore, it can be seen that the placement of the different particle sizes also has a large effect on the charging time. Since Configuration III has small particles closer to the inlet, a lot of energy can be absorbed from the HTF by the small particles close to the inlet. At a later time, Configuration IV overtakes Configuration III in terms of stored energy. This is due to the fact that smaller particles close to the inlet are already fully molten and thus not absorbing more heat. Therefore, Configuration III only has large particles left to absorb heat. Configuration IV on the other hand has its energy storage more spread out over the small and large particles. Therefore, all particles can still absorb heat at this stage, and, since small particles are more effective than large ones, the energy stored increases more rapidly. This effect is shown more clearly in Section 4.2.
The outlet temperature of the system is an indicator for the effectiveness of heat absorption by the bed. The temperature of the outlet is weighted by the mass flow through the n outlet CFD nodes over which the temperature is averaged. Since the HTF is considered to have a constant density, this leads to: The presented results in Figure 8 show that the temperature at the outlet for Configuration I is the highest. This is shown in Figures 6 and 7; since the HTF in Configuration I has transferred the least amount of energy to the particles, it leads to the most inefficient heat transfer and the highest outlet temperature. After about 3.5 h, the beds with smaller particles start overtaking Configuration I because, at that point, the small particles are already at the same temperature as the HTF. This means that less PCM is available to absorb heat from the HTF. After this point, Configuration I absorbs more energy since it is still in a less charged state (see Figure 7).

Effect of Placement Location of Small Particles
Another important question in designing a bed is on distribution and placement of particles with different sizes. To evaluate this effect, two cases are considered. The first one is Configuration IVa with large particles near the inlet and small particles near the outlet. In the other case (i.e., Configuration IVb), small particles are close to the inlet while large particles are near the outlet. These configurations are then studied during both charging and discharging.
As for the other cases, the average particles temperature is measured (see Figure 9). During both charging and discharging, the temperature lines of the two cases cross each other a couple of times.
Here, the charging and discharging cases mirror each other, although the discharging process is faster than the charging one. This is due to the melting temperature of the PCM being closer to the hot HTF temperature of 343.15 K as opposed to the cool HTF, being 293.15 K. Therefore, the temperature flatting of the PCM due to phase change creates a larger temperature difference with the HTF and, hence, a greater heat transfer during discharge. This process is also displayed in Figure 10. At the start of the charging, the temperature of Configuration IVb rises more quickly than in Configuration IVa, as shown in Figure 10a,b. This is due to the small particles being closer to the hot side of HTF, thus more heat transfer occurs. After a while, the lower rate of heat transfer to the large particles in Configuration IVa causes the HTF downstream to be hotter compared to Configuration IVb. Therefore, the small particles downstream in Configuration IVa which have already been heated up significantly lead to a rise in the average temperature.   Since latent heat storage has a more significant effect on the amount of stored energy (see Figure 11), more energy can be stored in the case of Configuration IVb for a duration of 2 h when compared to Configuration IVa, even though the average temperature of particles are lower in this case. This is shown in Figure 12a,b: it can be noted that, where particles in Configuration IVa are just undergoing phase transition, some particles in Configuration IVb are already fully molten, storing a lot of energy in the form of latent heat.
When the small particles in Configuration IVb move further along in the phase transition and are fully molten, the temperature starts to rise again and the average temperature of Configuration IVb passes Configuration IVa again, as can be seen in Figure 10c,d. This can also be seen by looking at the melt fraction of the particles in Figure 12c,d. Here, all the small particles in Configuration IVb, which is half of the total mass of the bed, have already undergone full phase transition and are heating up again. On the other hand, only a few large particles in Configuration IVa have passed phase change and a portion of the smaller ones. Hence, even though the temperature of Configuration IVb is higher, Configuration IVa can store more energy at this time. Eventually, the average temperature of Configuration IVa passes Configuration IVb while Configuration IVa remains the case which stores the most energy all the way to the end. This can be seen in Figure 12e,f. Here, most of the particles in Configuration IVa have fully undergone phase change so the temperature is rising again and the bed is almost fully charged. The bed of Configuration IVb still has more energy to store as latent heat. In the end, the bed of Configuration IVa achieves a fully charged state after 4 h 38 min while the bed of Configuration IVb reaches this state after 5 h 50 min. In the case of cooling, Configuration IVa is fully discharged after 2 h 13 min and Configuration IVb after 1 h 57 min. One can see that, during charging, Configuration IVb charges about 1 h more slowly than IVa, which is a decrease in charging time of 20%. During cooling, it discharges about 15 min more quickly, which is a decrease in charging time of 13%. Due to the larger temperature difference during phase change, the effect of the configuration is less present. Figure 13 shows that the average outlet temperature follows the same trends as the average particle temperature (see Figure 9). Initially, the energy transfer to Configuration IVb is higher due to the high energy transfer to the smaller particles but is surpassed later on by Configuration IVa when the smaller particles downstream start phase changing, until Configuration IVa is fully charged.

Conclusions
To investigate the effect of bed structure on the charging and discharging of packed bed of PCM, the open-source version of CFDEM R coupling was further developed to calculate the phase transition of PCM particles. This model was applied to multiple different configurations to study the effect of using different sizes of particles. The placement of smaller size particles with respect to the flow of HTF was considered. The rates of heat storage for all configurations were evaluated, revealing the following conclusions: • Smaller particles increase the rate of heat transfer between PCM and HTF; in the configurations considered, a decrease of 26% from Configuration I with the biggest particles and Configuration IV with the most small particles was found. • Using smaller particles intermixed with larger ones cause the average temperature of the bed not to show the flatting effect of phase change profoundly, as the smaller particles go through phase change more quickly.

•
Adding smaller particles downstream of the inlet has more effect on the rate of charging than placing them upstream. Placing the smaller particles downstream in the configurations discussed here, a decrease of 20% in charging timecan be achieved and 13% in discharging time.

•
Addition of smaller particles upstream can lead to a greater charge rate and total charge if heat flow is cut off before a certain point.
To further improve the predictions, It is interesting to add a model of a solar boiler before the inlet of the model and thus introduce a more realistic heat supply such as an operational TES system would experience. This could give a greater insight into the influence of small particle placement within an unstructured bed.

Abbreviations
The following abbreviations are used in this manuscript: