Microsoft Word - PRES22_0064.docx DOI: 10.3303/CET2294076 Paper Received: 08 April 2022; Revised: 21 April 2022; Accepted: 30 April 2022 Please cite this article as: Zálešák M., Charvát P., Klimeš L., Duda J., 2022, Assessment of Modelling Approaches for Partial Phase Charges of PCMs, Chemical Engineering Transactions, 94, 457-462 DOI:10.3303/CET2294076 A publication of CHEMICAL ENGINEERING TRANSACTIONS VOL. 94, 2022 The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Petar S. Varbanov, Yee Van Fan, Jiří J. Klemeš, Sandro Nižetić Copyright © 2022, AIDIC Servizi S.r.l. ISBN 978-88-95608-93-8; ISSN 2283-9216 Assessment of Modelling Approaches for Partial Phase Changes of PCMs Martin Zálešáka,*, Pavel Charváta, Lubomír Klimešb, Jiří Dudac a Department of Thermodynamics and Environmental Engineering, Brno University of Technology, Technická 2896/2, 61669 Brno, Czech Republic b Sustainable Process Integration Laboratory—SPIL, NETME Centre, Brno University of Technology, Technická 2896/2, 61669 Brno, Czech Republic c NUKLEON s.r.o., Ptáčov 40, 674 01 Třebíč, Czech Republic 162098@vutbr.cz Energy storage is going to play an important role in the transition to carbon-neutral energy production. Heat is the form of energy that is relatively easy to store, even over long periods of time. Despite their relatively slow take off as a thermal energy storage medium, the phase change materials (PCMs) can play an important role in the future of thermal energy storage (TES). Some PCMs exhibit behaviours that are difficult to model accurately (such as supercooling or phase change hysteresis). Computer modelling of the PCM behaviour is further complicated by the fact that the PCMs sometimes do not undergo complete phase transitions in real-life operation of TES. Three approaches to modelling of partial phase changes of a PCM, obtained from the literature, were assessed on a simple test case. The test case was a wall with a PCM layer, which was designed as transient with the adiabatic boundary condition on one side of the wall and the sine wave heat flux on the other side. The phase change hysteresis shift of 5 °C was investigated. The results obtained with the partial phase change modelling approaches varied on average by 2.8 °C and 1.84 °C for the ‘curve switch’ model and the ‘curve scale’ model (when compared to the single curve model). The increase in thermal conductivity (from 0.2 to 0.3 W m-1 K-1) led to 34.4 % reduction in size of the phase change interruption point interval. 1. Introduction Since most renewable energy sources are intermittent in their nature (reliant on weather conditions), there is an increasing need for means of energy storage, which would mitigate the energy supply/demand mismatch (Elfeky et al., 2021). Phase change materials (PCMs) are one of promising technologies for thermal energy storage (TES), mainly due to their ability to absorb and release a large amount of heat during the phase transitions (Serikawa et al., 2021). Some PCMs, however, exhibit behaviour such as phase change hysteresis or supercooling, which are difficult to model accurately. The computer modelling of PCMs gets complicated by the fact, that PCMs not always undergo complete phase transitions in TES systems. In case that phase transition is interrupted inside the ‘mushy’ region (region in-between the solid and liquid state) it is often referred to as a partial phase transition. Such cases are not uncommon, since most of the PCM-based TES systems are designed for operation around the phase change temperature of a PCM (Zastawna-Rumin et al., 2020). In recent years, a number of researchers addressed partial phase transitions in their studies (Klimeš et al., 2020). Thermal behaviour of a PCM is usually described by the relationships between temperature and enthalpy or temperature and effective heat capacity as is shown in Figure 2. Most modelling techniques confine the thermal behaviour of a PCM between the heating and cooling curves. The most straightforward approach is referred to as the ‘curve track’ model (Thonon et al., 2021), in which thermal behaviour (phase transition) is modelled by following either the heating or cooling curve. The model ‘tracks’ the heating or cooling curve until the material is considered fully molten (liquid) or fully congealed (solid) (Chandrasekharan et al., 2013). Close to the complete phase transition points, the transition between the curves can be made since there is little to no difference in enthalpy between the heating and cooling curves (Zalesak et al., 2021). The approach suggested by (Bony and Citherlet, 2007) utilises a linear transition with the slope derived from specific heat capacity 457 (see Figure 2, 2->3 transition), which is also very often referred to as the ‘curve switch’ model. However, recent studies suggest that the transition somewhere in-between the ‘curve track’ model and the ‘curve switch’ model is very promising (Barz, 2021). This modelling approach is called the ‘curve scale’ model, where the partial phase transitions follow the scaled cooling or heating curve (after interrupted heating or cooling process) (see Figure 3, the ‘curve scale’ model). The main objective of the present paper is to compare these three modelling approaches and to investigate the influence of thermal conductivity. 2. Methods The test case was a very simple one, it was a vertical layer (a wall) made of a PCM (Figure 1). The thickness of the wall was considered d = 100 mm and the PCM was paraffin-based Rubitherm RT35 HC. The enhancement of a building structure with a PCM for TES is a very common PCM application, allowing for shifting of the peak cooling load and improving indoor thermal comfort (Serikawa et al., 2021). The main advantage of building components enhanced with PCMs is their high heat capacity for working temperatures close to the phase change temperature. The latent heat storage therefore adds to the overall energy density. Figure 1: The PCM wall schematics 2.1 Numerical model The model in this study is constructed as a 1-D heat transfer problem (Figure 1). Heat conduction is considered only in the direction of PCM’s layer thickness and the governing equation is defined as 𝜌ceff ∂T ∂𝜏 = ∂ ∂𝑥 𝑘 ∂T ∂𝑥 (1) where T is temperature, 𝜏 stands for time, 𝜌 is density, 𝑘 represents thermal conductivity, ceff stands for effective heat capacity and x is the respective axis of the Cartesian coordinate system. It has been shown in several experimental studies (Serikawa et al., 2021), employing differential scanning calorimetry (DSC) method, that the effective heat capacity curve has an asymmetrical shape. The effective heat capacity in the present study was parametrised as ceff = c + c𝐿𝑓exp ― (T ― Tppc) 2 𝜎 (2) where c represents the specific heat capacity of solid and liquid phase, c𝐿𝑓 stands for the height coefficient (defining the amount of latent heat), Tppc is the peak phase change temperature and 𝜎 is the sharpness coefficient of ceff function defined as 𝜎 = 𝜎𝑠, 𝑇 ≤ 𝑇ppc 𝜎l, 𝑇 > 𝑇ppc (3) where 𝜎𝑠 and 𝜎𝑙 are representing the sharpness of the effective heat capacity curve on the heating and cooling sides. 458 The energy balance approach was employed for the solution of Eq(1), using the 1-D forward time-explicit numerical scheme. The phase change was modelled with the use of the enthalpy method, mainly due to its better numerical stability. The enthalpy method is also much more suitable for partial (incomplete) phase changes, since the phase transition is represented by a continuous function; contrary to the effective heat capacity method (as illustrated in Figure 2). The volume enthalpy is defined as 𝐻(𝑇) = 𝑇 𝑇ref 𝜌c ― 𝜌𝐿f ∂𝑓s ∂𝜃 d𝜃 (4) where 𝐿𝑓 is the latent heat of phase change and 𝑓𝑠 represents the melt fraction. The calculation process turns into the 2-step procedure – in the first step, the enthalpy of (𝑝 + 1)-th time discretization step (𝐻p+1) is calculated. According to knowledge of the enthalpy-temperature dependence, the temperature 𝑇p+1 is then determined from the known enthalpy-temperature dependence. The entire process can be summarised as 𝑇p→ 𝐻p+1→ 𝑇p+1. Figure 2: Effective heat capacity and enthalpy method considering the “curve switch” model The boundary conditions on both sides of the wall are of Neumann’s type. On the left (Figure 1) the adiabatic boundary condition was prescribed (𝑞𝐿 = 0). On the right (Figure 1), the wall was exposed to the cyclical heat flux 𝑞𝑅 formulated as 𝑞𝑅 = 𝐴0sin 2𝜋𝑖 𝑝max ,∀𝑖 ∈ 〈0, 𝑝max〉 (5) where 𝑖 is the time instance and 𝑝max is the total number of iterations. 2.2 The studied scenario Thermophysical properties of Rubitherm RT35 HC, together with other considered parameters, are summarised in Table 1. The total duration of the simulation was ∆𝑡max = ∆𝑡 ∙ 𝑝max = 0.1 ∙ 30,000 = 3,000 s. In time 𝑡 = 0 s, the PCM is fully in a solid state with the initial condition 𝑇p =0 𝑗 = 15, ∀𝑗 ∈ 〈1,𝑁𝑥〉. The scenario is designed for the investigation of partial melting-solidification transitions. The amplitude of the sine boundary condition 𝐴0 = 16 W m2 (see Equation 5) was chosen in a way, which ensures that the PCM in all nodes of the computational domain undergoes melting-solidification transitions inside the phase change temperature interval (partial phase transitions). Heat losses into the surroundings were not considered. The solidification curve was shifted by ∆𝑇𝐻 = 5 °C to lower temperature as illustrated in Figure 2. As for the spatial discretization of the wall, in total 𝑁𝑥 = 20 nodes were considered with ascending indexes starting from the surface with the prescribed sine wave boundary condition (see Figure 1). 459 Table 1: Properties of Rubitherm RT 35 HC and other model parameters Material property Symbol Unit Value Phase change temperature Tppc [°C] 35 Phase change hysteresis shift ∆𝑇𝐻 [°C] 5 Heat storage capacity* 𝐿𝑓 [kJ/kg] 210 Specific heat capacity c [J/(kg K)] 2,000 Time discretization step ∆𝑡 [s] 0.1 Density 𝜌 [kg/m3] 800 Thermal conductivity 𝑘 [W m-1 K-1] 0.2 Number of nodes 𝑁𝑥 [-] 20 Wall thickness d [m] 0.1 Heating/cooling switch threshold THR [-] 50 *Excluding the sensible heat between 27 - 42 °C 3. Results Based on the literature review, three most commonly used modelling approaches to phase change hysteresis of PCMs were chosen for testing in the present simulation study; the ‘curve track’ model (as well referred to as no transition model or single curve model), the ‘curve switch’ model (implementation proposed by Bony and Citherlet, (2007)) and the ‘curve scale’ model. The sine wave heat flux boundary condition gradually heated up the whole volume of a PCM into the phase change interval. As soon as there were 𝑁𝑖𝑡ℎ ≥ THR = 50 consecutive temperature gradients (in a i-th node) monitored, the heating-to-cooling transition was performed following one of the implemented approaches. As can be seen in Figure 3, the separation point of the node 3 did occur at about 21 min with temperature of 34.8 °C, which is well within the phase transition interval and it is marked as the transition point in this specific scenario and node. Same observations can be made about other nodes as they all ended up in the phase change temperature interval of the PCM (see Figure 3, nodes 8, 13, 18). Thermal behaviour of the PCM after the transition point was influenced by the particular modelling approach to phase change hysteresis and the discrepancies (the temperature differences between the specific partial phase change modelling approach and the no transition model) were on average 2.8 °C and 1.84 °C for the ‘curve switch’ model and the ‘curve scale’ model. The ‘curve track’ model (no transition) seems to be overestimating the temperature evolution in comparison to the other two models (see Figure 3), with the ‘curve scale’ model representing the approach which is in-between the no transition (the ‘curve track’ model) and the linear transition between the curves (the ‘curve switch’ model). This behaviour is in agreement with the previous research presented by (Delcroix et al., 2017), suggesting that the ‘curve scale’ modelling approach could be potentially very promising. Figure 3: Temperature evolution of the implemented phase change hysteresis modelling approaches 460 Material properties such as thermal conductivity, specific heat capacity or enthalpy of fusion affect the thermal behaviour during the phase change hysteresis. The results have shown that with the increasing thermal conductivity of a PCM the interval of the heating-to-cooling transition points is narrower, since the overall temperature gradient inside the PCM wall decreases (as is illustrated in Figure 4). On the other hand, in case of a PCM with low thermal conductivity, fewer nodes would be in the mushy zone, under the same boundary conditions, as the temperature gradient in the PCM would be steeper. At the heated surface, PCM could be fully molten while PCM at the insulated surface could still be in the solid state (depending on heat flux rate). The increase in thermal conductivity from 0.2 to 0.3 W m-1 K-1 resulted in 34.4 % decrease in size of the interval of phase change interruption points (from [33.51, 35.34 ] °C to [33.69, 34.89] °C) as is shown in Figure 4. Figure 4: The effect of thermal conductivity on phase change hysteresis 4. Conclusions The present study aimed at the assessment of recently proposed modelling approaches to the phase change hysteresis phenomenon and the sensitivity analysis with regard to thermal conductivity of PCMs. The 1-D numerical model of heat transfer in PCM was developed using the enthalpy method and it was coupled with three different partial phase change modelling approaches. There was a considerable overestimation of the temperature when the ‘curve track’ model was compared with the ‘curve switch’ and the ‘curve scale’ models, suggesting that neglecting of the phase change hysteresis (with phase change temperature shift ∆𝑇𝐻 = 5 °C) can result in the average temperature discrepancies of 2.8 °C and 1.84 °C. With the increasing thermal conductivity, heating-cooling transition points converge as the overall temperature gradient in the wall gets reduced. The results have shown that increase in thermal conductivity (from 0.2 to 0.3 W m-1 K-1) led to 34.4 % reduction of the size of the phase change interruption point interval. On the other hand, decreasing thermal conductivity leads to complete phase transitions at the heated surface and the insulated surface might not even fully undergo the phase change. Therefore, the phase change transition point interval is getting wider. The present study provides the foundation for inverse identification of thermal behaviour of PCMs during the partial phase changes, as it is difficult to investigate thermal behaviour of large volumes of PCMs with the conventional experimental methods such as DSC or T-history (which are designed for small PCM samples). The future work will consider the potential effects of other material properties such as specific heat capacity or latent heat of fusion on the thermal behaviour of PCMs during partial phase changes. The experimental validation of the partial phase change modelling approaches is also planned for the future studies. As the melting front becomes much more complex with the higher dimensionality of the numerical model (considering 2-D and 3-D models), other aspects such as natural convection in the liquid PCM or the volumetric change (caused by the change of density between the solid and liquid states of PCMs) have to be accounted for. Nomenclature 461 𝐴0 – amplitude of the sine wave boundary condition, W m2 𝑐 – specific heat capacity, J/(kg K) c𝐿𝑓 – hight coefficient of the ceff, J/(kg K) ceff – effective heat capacity, J/(kg K) d – PCM wall thickness, m 𝑓𝑠 – melt fraction, - 𝐻(𝑇) – volume enthalpy, J 𝑘 – thermal conductivity, W/(m K) 𝐿𝑓 – enthalpy of fusion, J/kg 𝑁𝑥 – number of nodes, - 𝑝max – total number of iterations, - T – temperature of PCM, K Tppc – peak temperature of phase change, K THR – heating/cooling switch threshold x, y – cartesian coordinates ∆𝑡 – time discretization step, s ∆𝑡max – total duration of simulation, s ∆𝑇𝐻 – phase change hysteresis shift, K 𝜌 - density of PCM, kg/m3 𝜎 – sharpness coefficient of the ceff, K2 Acknowledgements This work was supported by the project Computer Simulations for Effective Low-Emission Energy Engineering, reg. no. CZ.02.1.01/0.0/0.0/16_026/0008392 funded by Operational Programme Research, Development and Education, by the project Adaptive Soft Computing Framework for Inverse Heat Transfer Problems with Phase Change, reg.no. 22-31173S funded by the Czech Science Foundation, and by the internal research project of Brno University of Technology, reg. no. FSI-S-20-6295. References Barz T., 2021. Paraffins as phase change material in a compact plate-fin heat exchanger - Part II: Validation of the ‘‘curve scale’’ hysteresis model for incomplete phase transitions. Journal of Energy Storage, 34, 102164. Bony J., Citherlet S., 2007. Numerical model and experimental validation of heat storage with phase change materials. Energy and Buildings, 39, 1065-1072. Chandrasekharan R., Lee E.S., Fisher D.E., Deokar P.S., 2013. An enhanced simulation model for building envelopes with phase change materials. ASHRAE TRANS, 119, 1-10. Delcroix B., Kummert M., Daoud A., 2017. Development and numerical validation of a new model for walls with phase change materials implemented in TRNSYS. J Build Perform Simul, 10, 422–37. Elfeky K.E., Mohammed A.G., Wang Q., 2021. Numerical Study to Investigate the Thickness of the PCM Layer for the Three Layers Tank for the CSP Plants, Chemical Engineering Transactions, 88, 85-90. Klimes, L., Charvat, P., Joybari, MM., Zalesak, M., Haghighat, F., Panchabikesan, K., El Mankibi, M., Yuan, YP., 2020. Computer modelling and experimental investigation of phase change hysteresis of PCMs: The state-of-the-art review. Applied Energy, 263, 114572. Serikawa M., Mabuchi K., Satoh M., Nozue Y., Hayashi Y., Yokoyama M., 2021. Measurement of full-scale phase change material products considering hysteresis. Applied Thermal Engineering, 192, 116895. Thonon, M., Fraisse, G., Zalewski, L., Pailha, M., 2021. Analytical modelling of PCM supercooling including recalescence for complete and partial heating/cooling cycles. Applied Thermal Engineering, 190, 116751. Zalesak M., Charvat P., Klimes L., 2021. Identification of the effective heat capacity-temperature relationship and the phase change hysteresis in PCMs by means of an inverse heat transfer problem solved with metaheuristic methods. Applied Thermal Engineering, 197, 117392. Zastawna-Rumin, A., Kisilewicz, T., Berardi, U., 2020. Novel Simulation Algorithm for Modeling the Hysteresis of Phase Change Materials. Energies, 13, 1200. 462