- Email: [email protected]
Simultaneous charging and discharging of phase change materials: Development of correlation for liquid fraction
Solar Energy 188 (2019) 788–798
Contents lists available at ScienceDirect
Solar Energy journal homepage: www.elsevier.com/locate/solener
Simultaneous charging and discharging of phase change materials: Development of correlation for liquid fraction
T
⁎
Mahmood Mastani Joybaria, Fariborz Haghighata, , Saeid Seddeghb, Yanping Yuanc a
Energy and Environment Group, Department of Building, Civil and Environmental Engineering, Concordia University, Montreal H3G 1M8, Canada School of Engineering, College of Science and Engineering, University of Tasmania, Hobart, TAS 7001, Australia c School of Mechanical Engineering, Southwest Jiaotong University, Chengdu 610031, China b
A R T I C LE I N FO
A B S T R A C T
Keywords: Front tracking Phase change material Natural convection Triplex tube heat exchanger Simultaneous charging and discharging
In this study, a method to track the melting front of phase change materials (PCMs) is developed under simultaneous charging and discharging (SCD) inside horizontal triplex tube heat exchangers. For SCD, the heat transfer mode of internal heating and external cooling is investigated. According to the heat transfer mechanism within the PCMs, two models of pure conduction (PC) and combined conduction and natural convection (CCNC) are developed and solved using ANSYS Fluent (version 17.2). Generally, during charging (i.e. melting), the buoyancy forces induce the liquid PCM to move upwards; therefore, the upper section of the storage is affected. However, this motion is limited during SCD due to the simultaneous cooling from the outer tube. Therefore, a front tracking method is developed assuming that natural convection only contributes melting to the upper section of the storage unit. On the other hand, it is assumed that the melting front for the lower section behaves similar to the PC model. In this study, three PCMs are used in three geometries to develop the correlations using two dimensionless numbers; i.e. liquid fraction of the PC model and the storage radius ratio. Another PCM is then used to verify the method. The correlations can provide results in the range of ± 15% discrepancy.
1. Introduction Latent heat storage using phase change materials (PCMs) benefits from high storage capacity and can adjust the time mismatch between energy supply and demand (Mirzaei and Haghighat, 2012). Therefore, it found several applications, including building envelopes (Liu et al., 2018), ventilation systems (El-Sawi et al., 2014), refrigeration systems (Joybari et al., 2015), net zero energy buildings (Stritih et al., 2018), hot water tanks (Najafian et al., 2015), etc. Prior to designing any latent heat storage system using PCMs, the heat transfer mechanisms should be understood. In the literature, heat transfer within PCMs can be modeled using (i) the pure conduction (PC) model or (ii) combined conduction and natural convection (CCNC) model. The PC model was used in earlier studies to investigate the heat transfer within PCMs (Seddegh et al., 2015). Since gravity has no effect on conduction, the melting front would be circular (or cylindrical in 3D) (Tay et al., 2012a; Joybari and Haghighat, 2016). In reality, however, as the PCM melts, buoyancy forces due to the change in PCM density create an upward motion (Seddegh et al., 2016). This
phenomenon is known as natural convection, which has been observed in several experimental studies (Al-Abidi et al., 2013b; Joybari et al., 2019). The two main design tool categories for heat exchangers are (i) simplified methods, and (ii) complex methods including CFD simulations. The first category includes the well-known approaches such as effectiveness-number of transfer units (ε-NTU), normally oversimplifying the heat transfer process within PCMs. For instance, neglecting natural convection results in effectiveness underestimation during the phase change process (Tay et al., 2012a,b,c; Bruno, 2008). In the second category, CFD simulations are more complicated and normally require long computational time. Nevertheless, some tools can be developed, whose simplicity and accuracy lie between the two mentioned categories. To enhance the accuracy of the ε-NTU method, a method called effective thermal conductivity was introduced which accounts for the natural convection during the melting process. In this method, the thermal conductivity values are artificially increased to mimic the experimental results. Comparison of ε-NTU method with experimental
Abbreviations: CFD, computational fluid dynamics; CCNC, combined conduction and natural convection; HTF, heat transfer fluid; NTU, number of transfer units; PC, pure conduction; PCM, phase change material; STHX, shell and tune heat exchanger; TTHX, triplex tube heat exchanger ⁎ Corresponding author. E-mail address: [email protected] (F. Haghighat). https://doi.org/10.1016/j.solener.2019.06.051 Received 6 March 2019; Received in revised form 28 May 2019; Accepted 20 June 2019 0038-092X/ © 2019 International Solar Energy Society. Published by Elsevier Ltd. All rights reserved.
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
density (kg m−3)
Nomenclature
ρ
a b c C∗ Fo h∗ P* Pr r R Ra Ste V
Subscripts
constant constant constant dimensionless Darcy coefficient (-) Fourier number (-) dimensionless enthalpy (-) dimensionless pressure (-) Prandtl number (-) radius (m) storage radius ratio (-) Rayleigh number (-) Stefan number (-) dimensionless velocity (-)
Superscripts final L section modified U section dimensionless variable
f L mod U ∗
Greek symbols
γ
combined conduction and natural convection model inner outer pure conduction model reference
CCNC i o PC ref
liquid fraction (-)
In this method, both phases (i.e. solid and liquid) are modeled using the same enthalpy-based energy equation (Voller and Prakash, 1987). Therefore, to track the melting front using the enthalpy method, a twostep process should be followed: (i) calculation of the enthalpy values, and (ii) determination of the melting front location according to the temperature values (Li et al., 2003). Nevertheless, despite the progress in computational power, the numerical analysis of PCMs is still complicated. It also requires code development expertise or accessibility to a commercial software for simulation. Therefore, to characterize the PCM phase change process, several correlations have been reported in the literature (see Table 1). As the table indicates, the older correlations were developed for Fourier number (Fo) to obtain the phase change duration. However, more recent studies focused on liquid fraction (γ). Note that natural convection was not considered in most of the studies. Overall, although correlations have been developed for PC and CCNC models, to the best of the authors’ knowledge, no correlation has been developed linking these two models. Another inference from Table 1 is that PCMs have normally been investigated under one-time (or consecutive processes of) charging and/or discharging. In other words, an initially solidified (or melted) PCM was considered to undergo melting (or solidification) process. This ideal assumption is not in line with real implementations where the PCM might undergo simultaneous charging and discharging (SCD). In such cases, a PCM is being melted from one side, while it is
data revealed that the method would have acceptable results only if it is enhanced by the effective thermal conductivity method (Tay et al., 2012c). Nevertheless, there are some drawbacks for the effective thermal conductivity method: (i) it requires experimental data; therefore, an experimental setup should be configured and manufactured, (ii) it is time-consuming as the best value is obtained by trial and error while comparing to experimental data, and (iii) despite all these complexities, since the effective thermal conductivity method is based on the conduction concept, no information can be obtained regarding the melting front. Knowing the location of melting front can not only indicate the liquid fraction value, but also capture partial melting. However, obtaining liquid fraction in experimental studies is very complicated. During earlier experimental studies the phase change process was interrupted after some time in order to extract the remaining solid part (Sparrow and Broadbent, 1982). In recent years, transparent tubes were utilized which can provide direct visual observation (Liu and Groulx, 2014). Moreover, digital high resolution photography (Yang et al., 2016) as well as image processing (Jones et al., 2006) has greatly facilitated the liquid fraction calculation process. Nevertheless, it is preferable to conduct numerical analyses in order to obtain liquid fraction. Melting front tracking used to be a great challenge in early numerical studies since a moving-boundary problem (due to the change of melting front location) used to be solved (Klimeš et al., 2018). However, such a complexity can be avoided by using the enthalpy method.
Table 1 Hierarchical summary of the existing correlations characterizing the phase change process. Refs.
Type of study
Storage type(s)
Orientation
(Riley et al., 1974)
Analytical
CylindricalSpherical
Fo = a +
b Ste
(Solomon, 1979)
Analytical
Fo = a +
b Ste
(Hasan, 1994b)
Experimental
Slab CylindricalSpherical Cylindrical
Fo = a +
b Ste
(Hasan, 1994a)
Experimental
Cylindrical
Fo = a +
b Ste
(Ismail and Maria das Graças, 2003)
Numerical
Cylindrical
Horizontal Vertical Horizontal Vertical Horizontal
(Katsman et al., 2007)
ExperimentalNumerical
STHX
Vertical
γ = (FoSte )aRaDb exp
γ = aFoSte bRac Fo = a + b Re + cSte
(Fan et al., 2014)
Experimental
Cylindrical
Vertical
(Rathod and Jyotirmay, 2015)
Experimental
STHX
Vertical
789
Correlation
Comment Developed for discharging (no significant natural convection occurred). Developed for both charging and discharging.
Y = a (X )2 + b (X ) + c
( ) c RaH
Developed for both charging and discharging. Natural convection was neglected despite reporting its effect. Developed for both charging and discharging. Natural convection was neglected despite reporting its effect. Developed for charging. The correlations include complete melting duration and melted volume Stefan (Ste) or Rayleigh (Ra) number is represented by X. Developed for charging. Developed for charging. Developed for charging.
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
• No correlation has been developed for PCMs linking the CCNC and
simultaneously being solidified from the other. This can happen when a PCM is unpredictably and/or intermittently charged and/or discharged. As an example, when PCM is used for solar domestic hot water, occupants might demand water consumption (i.e. discharging) during sunshine (i.e. charging) hours. According to the literature, despite several studies which investigated the charging, discharging or consecutive charging and discharging process, the behavior of PCMs when undergoing SCD has not been thoroughly studied. For instance, Liu et al. (2006) experimentally investigated the SCD process in a heat pipe heat exchanger. A paraffin PCM was used where its high thermal resistance was found to inhibit efficient system performance. However, that was not the case for a metal high-temperature PCM in another study (Sharifi et al., 2015). Murray and Groulx (2014) conducted experimental investigations for SCD of a PCM for domestic hot water applications. Charging was continuously imposed, while the discharging was intermittent. Again, the solid part of the PCM controlled the heat exchange process due to its high thermal resistance. Omojaro and Breitkopf (2014) numerically and experimentally investigated SCD in a horizontal shell and tube heat exchanger (STHX) where the PCM was located in the shell. For charging and discharging, a heater and a fan were used, respectively. For numerical investigation, a simple symmetric one-dimensional model was developed for temperature as well as liquid fraction. Recently, the same authors used the existing experimental setup to investigate the system performance under different convection mechanisms from the shell (i.e. natural convection and forced convection) (Omojaro and Breitkopf, 2017). It was reported that forced convection resulted in lower melting rate. Besides, lower air velocities caused shorter transition time within the PCM from conduction to convection. Thus, the low thermal conductivity of the PCM deteriorated the system performance. Certain configurations have the capability to undergo SCD. Triplex tube heat exchangers (TTHXs) consist of three concentric tubes. For heat storage, the PCM fills the middle space, while the HTFs flow through the inner and outer tubes. TTHXs also possess some characteristics, suitable for solar domestic hot water applications: (i) possibility to impose SCD, (ii) possessing separate tubes for the HTF and water to flow which can satisfy sanitary requirements, and (iii) enhanced natural convection, improving heat transfer within the solid PCM. Several studies investigated the phase change process in PCMequipped TTHXs under one-time solidification (Al-Abidi et al., 2013e; Abdulateef et al., 2018; Mahdi and Nsofor, 2018), one-time melting (AlAbidi et al., 2013b, 2013c, 2013d; Mat et al., 2013), or consecutive melting and solidification (Al-Abidi et al., 2014; Jian-you, 2008; Long and Zhu, 2008). In addition, SCD in TTHXs has recently gained attention. For instance, Mahdi et al. (2019) recently optimized the fin configuration in TTHXs under SCD and reported that it outperformed addition of nanoparticles (i.e. a PCM thermal conductivity enhancement technique). In our earlier study, a horizontal TTHX was numerically investigated under SCD (Joybari et al., 2017c). For internal heating/ external cooling, it was found that the PCM can be partially melted in the upper half of the system. Later, fin configurations were recommended under SCD as well as non-SCD applications for TTHXs in accordance with the behavior of natural convection (Joybari et al., 2017a). Reviewing the available literature demonstrates the following limitations:
• • • •
PC models under SCD.
Therefore, in this study, a simple method of front tracking is developed for SCD, considering the natural convection behavior. In our earlier study (Joybari et al., 2017b), a method was developed for melting process only. Therefore, this study extends its application to SCD providing important information regarding the melting front location and partial melting. 2. Materials and methods The considered storage geometrical properties as well as PCM thermophysical properties are described in this section, which are the same as our earlier study for melting (Joybari et al., 2017b). A summary of the modeling approach is also presented. 2.1. Setup characteristics The US standard for sizing of copper tubes was used which recommends tubes of “Type M” for water heating applications (Handbook et al., 2011). A horizontal TTHX was investigated in which the PCM filled the middle tube (see Fig. 1). Three geometries (see Table 2) were used for data generation (to be used for correlation development). Note that the inner tube radius (ri) was altered, while the other tube size (ro) remained the same. 2.2. PCM thermophysical properties In this study, four commercially available paraffin PCMs (from Rubitherm Technologies GmbH) were investigated (see Table 3). Note that during data generation process the highlighted PCM (i.e. RT35HC) was excluded so that it can be used later for correlation verification. 2.3. Governing equations This section discusses the assumptions as well as the modeling procedure adopted for this study. 2.3.1. Assumptions The following were considered for the phase change process of the PCMs:
• Assumed
laminar, incompressible Newtonian fluid flow for the
ri
Most earlier studies (especially older ones) neglected the natural convection during charging and/or discharging process, CFD simulations using the CCNC model are not only time consuming but also require high computational power, The phase change process under SCD has not been thoroughly studied in the literature and experimental or numerical investigations of PCM under SCD are limited, especially in TTHXs, There is no correlation in the literature for front tracking under SCD,
ro ș=0
Fig. 1. Schematic representation of the TTHXs. 790
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
Table 2 Geometrical properties of the middle space in the TTHXs (Joybari et al., 2017b; Handbook et al., 2011). Geometry
Geometry 1 Geometry 2 Geometry 3
Radius (mm)
Storage radius ratio (R)
Inner (ri)
Outer (ro)
26.98 39.69 52.38
62.32 62.32 62.32
2.31 1.57 1.19
melted PCM,
• Neglected viscous dissipation due to the absence of large velocity gradients, • Assumed no-slip and Dirichlet boundary conditions for the fluid flow and heat transfer, respectively, • Assumed linear variation of liquid fraction with temperature, • Assumed the thermophysical properties to be constant for each
Fig. 2. Model validation under the melting process.
solidification processes. The simulation results are compared with experimental data obtained from two earlier studies.
phase except density, being evaluated using Boussinesq approximation (considering natural convection).
2.3.4.1. Melting process. The experimental data of Mat et al. (2013) was used to validate the model under melting. The paraffin-based RT82 (Rubitherm Technologies GmbH) was used in that study as the PCM in a horizontal finned TTHX. To melt the PCM, hot water was used, flowing inside both inner and outer tubes. Experimental data (i.e. average temperature) is compared with numerical results in Fig. 2. Despite slight overestimation of the model, an acceptable agreement exists with the experimental data. Therefore, the model is considered validated under melting.
2.3.2. Modeling Our earlier studies (Joybari et al., 2017a,c) include the detailed modeling approach. For the sake of conciseness and generalization, this section presents the equations using dimensionless parameters: Continuity Equation:
→ ∂ρ + ∇ ·(ρV ) = 0 ∂Fo
(1)
Momentum Balance:
→ → →→ → → ∂V Ra ·Pr ∗ ∗ + ∇ ·( V V ) = −∇P ∗ + Pr ∇ ·(∇ V ) + (h − href ) − C ∗V ∂Fo Ste
2.3.4.2. Solidification process. Experimental data of Seddegh et al. (2017a) was used to validate the model under solidification. The paraffin-based RT60 (Rubitherm technologies GmbH) was used as the PCM in a vertical STHX. Water was used as the HTF in the tube to solidify the PCM. The temperature of an exact point (i.e. T8 of Cylinder D in Seddegh et al. (2017a) was considered for validation. Fig. 3 shows the results where the model followed the trend of experimental data. Note that the deviations are because any thermocouple records the temperature of its surrounding chunk of the PCM. In other words, they are not capable of measuring the point-specific temperature. Despite this, the numerical results reasonably agree with those obtained experimentally. Therefore, the model is considered validated under solidification as well.
(2)
Energy Balance: ∗ → ∂γ ⎡ ∂h + ∇ ·(→ + ∇ ·( V γ ) ⎤ = ∇ ·(∇h∗) V h∗) ⎤ + ⎡ ⎣ ∂Fo ⎦ ⎣ ∂Fo ⎦
(3)
2.3.3. Numerical approach Several commercial software exists for numerical investigation of the phase change process in PCMs. Since ANSYS Fluent (based on the enthalpy method) is most commonly used (Al-abidi et al., 2013a), it was adopted for this study. First, the 2D geometry was developed in DesignModeler. Thereafter, ANSYS Meshing was used to generate triangular meshes. Then, simulation setup was carried out in ANSYS Fluent. The pressure-based solver was utilized as it is compatible with solving phase change processes (ANSYS Fluent Theory Guide, 2013). For pressure-velocity coupling, the SIMPLE method was used. QUICK and PRESTO schemes were adopted to spatially discretize the momentum and pressure fields, respectively.
2.3.5. Analysis procedure In this study, the simulations were performed for each PCM (details in Table 3) in various geometries (details in Table 2). Consequently, Table 4 tabulates all resulting combinations for Stefan numbers as well as storage radius ratios. Note that in Table 4 three cases are designated as “Ver”. These indicate cases which are later used for correlation verification. Overall, 42 simulations were conducted as both the PC and CCNC models were simulated.
2.3.4. Model validation In this section, the model is validated separately for melting and
Table 3 The PCMs and their thermophysical properties (Joybari et al., 2017b; GmbH, 2017). Property
Dimension
RT31
RT35
RT35HC
RT44HC
Melting onset temperature Solidification onset temperature Solid density Liquid density Specific heat capacity Latent heat Thermal conductivity Thermal expansion coefficient Dynamic viscosity
K (°C) K (°C) kg/m3 kg/m3 J/kg K J/kg W/m K 1/K kg/m s
300.15 (27) 306.15 (33) 880 760 2000 170,000 0.2 0.00076 0.002508
302.15 (29) 309.15 (36) 860 770 2000 160,000 0.2 0.00076 0.002500
307.15 (34) 309.15 (36) 880 770 2000 240,000 0.2 0.00076 0.002700
314.15 (41) 317.15 (44) 800 700 2000 250,000 0.2 0.00076 0.003300
791
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
cooling from the outer tube, maintaining a solid PCM layer adjacent to the outer tube wall. 3.3. Comparison According to the results, compared to the CCNC model, the PC model benefits from:
• Easier analysis setup as well as shorter running time, • Easier meeting of convergence criteria for its energy equation (since continuity as well as momentum equations are not present), • Possibility of analytical solution (under certain circumstances), drastically reducing the computational time.
According to these, the present study is focused on developing a simple method, which can generate the CCNC model results using PC model data. The method can at the same time approximately locate the melting front, outperforming the effective thermal conductivity method.
Fig. 3. Model validation under the solidification process. Table 4 The considered cases indicating their Stefan number and storage radius ratio. Case
R
Ste
Case
R
Ste
Case
R
Ste
Case 1 Case 2 Case 3 Case 4 Case 5 Case 6 Ver 1
2.31 1.57 1.19 2.31 1.57 1.19 2.31
0.16 0.16 0.16 0.24 0.24 0.24 0.33
Case 7 Case 8 Case 9 Case 10 Case 11 Case 12 Ver 2
2.31 1.57 1.19 2.31 1.57 1.19 1.57
0.25 0.25 0.25 0.32 0.32 0.32 0.17
Case 13 Case 14 Case 15 Case 16 Case 17 Case 18 Ver 3
2.31 1.57 1.19 2.31 1.57 1.19 1.19
0.47 0.47 0.47 0.50 0.50 0.50 0.33
4. Performance analysis of the existing method In our earlier study (Joybari et al., 2017b), a method of front tracking was developed for melting in STHXs. The same geometry configurations as this study were used to develop the correlations. Fig. 6 shows the procedure used in that study for the method development. In this section, the possibility for application of the existing method under SCD is verified. According to Fig. 6, the existing method separated the melting front into two parts for each half of the storage based on the mechanism of natural convection. Since natural convection affects the upper half of horizontal storages, its effect was first assumed to be solely on the upper half. Consequently, the behavior of the lower half was considered identical to the PC model. However, after the upper half was melted completely, the existing method deviated the lower half by attributing the remainder of natural convection. To verify the performance of the existing method under SCD, three scenarios have been considered as shown in Table 7. The first two scenarios consider liquid fractions of 0.5 and 0.75 for the upper half, while Scenario 3 is for the steady state condition under SCD. In Fig. 7 the right and left sides of each contour show the results for the PC and CCNC models, respectively. According to the contours, since the upper half of the system never reaches complete melting under SCD (as shown in Fig. 5), the existing method (Joybari et al., 2017b) cannot accurately capture the mechanism of natural convection. Moreover, Fig. 8 presents the quantitative comparison between the existing method and the observed one. As the figure shows, the performance of the existing method improves towards the steady state condition (i.e. Scenario 3) mainly due to the deceleration of natural convection. Note that a thermal energy storage system should not always run under SCD. The reason is that under such conditions there would be no benefit for the PCM. In other words, the HTFs can have direct heat transfer. Consequently, an ideal method should be capable of providing more accurate results specially at the early stages of SCD (as the steady state condition may never occur). Therefore, the existing method may be applied under SCD only if steady state condition is to be reached; otherwise, it should not
Before conducting the analysis, the simulation independency in terms of grid size and timestep size were checked. The results for Case 4 are presented in this section as an example. To check grid size independency three resolutions (i.e. 2574, 4091 and 7930 elements) were considered. In addition, three timestep sizes (i.e. 0.05, 0.1 and 0.2 s) were also investigated. Tables 5 and 6 tabulate the liquid fraction results for these criteria under SCD for 80 min. According to the results, 0.1 s timestep with 4091 elements was found to be sufficient for the simulation. 3. Results and discussion Prior to method description, its underlying concept is elaborated by explaining and comparing the mechanism of PC and CCNC models. 3.1. Pure conduction model As its name implies, the PC model only considers conduction and neglects natural convection during melting. Although the PC model is not realistic, its simplicity and fast solution (under certain circumstances it could be analytically solved) make it favorable. Fig. 4 illustrates the melting front propagation for RT31 under SCD using the PC model. The liquid fraction contours show that the melting front is circular as it monotonically propagates in the radial direction. 3.2. Combined conduction and natural convection model According to experimental data (Seddegh et al., 2017b) when a PCM melts, buoyancy forces are formed due to the density change, creating natural convection. Fig. 5 illustrates how the melting front propagates for RT31 under SCD using the CCNC model. Comparing Fig. 5 (for the CCNC model) with that of the PC model (i.e. Fig. 4) it can be understood that the upward liquid PCM motion affects the upper section of the storage, enhancing heat transfer and subsequently shortening the phase change process. Another very important observation is that despite the upward motion, for none of the cases the upper half achieved complete melting (i.e. a liquid fraction of one). This is due to the simultaneous
Table 5 Verification of grid size independency (with 0.1 s timestep). Number of elements
2574 4091 7930
792
Simulation time (s) 300
600
1200
2400
4800
0.0238 0.0278 0.0312
0.0442 0.0540 0.0589
0.0918 0.1098 0.1153
0.1927 0.2283 0.2202
0.3632 0.3777 0.3354
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
part of the PCM (adjacent to the outer tube) will remain solid; thus, complete melting cannot be achieved. A closer look at Fig. 5 reveals that within the upper half of the system, the effect of natural convection on the top part is stronger than the bottom part. To further clarify this, an exemplary contour is shown in Fig. 9a. Despite simultaneous cooling, the upward melted PCM motion could melt those parts of the PCM denoted by U (i.e. the upper one third of the storage). The melted region in U is greater than the rest (denoted by L). This type of sectionalizing has been defined to meet method simplicity requirements upon scrutinized examination of liquid fraction contours under SCD. Consequently, natural convection is first attributed to propagate the location of melting front at U (see Fig. 9b) until its liquid fraction reaches a value of one. In other words, there would be two separate melting fronts where the front radius of U is larger than that of L. The developed correlation for U (i.e. f1) would map the results of the PC model to those of U. Once the PCM in U reaches complete melting (i.e. a liquid fraction value of one), the remainder of natural convection is designated to L. In other words, L starts to deviate from the PC model once U is completely melted. The developed correlation for L (i.e. f2) would map the results of the PC model to those of L. The phase change process would later stabilize once the steady state is reached under SCD. For the two melting fronts shown in Fig. 9, the overall liquid fraction is:
Table 6 Verification of timestep size independency (with 4091 elements). Timestep (s)
0.05 0.1 0.2
Simulation time (s) 300
600
1200
2400
4800
0.0276 0.0278 0.0279
0.0537 0.0540 0.0542
0.1083 0.1098 0.1101
0.2228 0.2283 0.2302
0.3667 0.3777 0.3836
be used, as it is not compatible with natural convection under SCD. 5. Front tracking method development for simultaneous charging and discharging Due to the improper performance of the existing method, a new method of front tracking together with its liquid fraction correlations are developed in this section specifically for SCD. Thereafter, a new PCM is used for verification of the correlations. 5.1. Description of the method When natural convection is absent (i.e. the PC model), the shape of melting front would be cylindrical (i.e. circular in 2D). Consequently, any fraction of the storage would have the same liquid fraction as the PC model. However, experimental data (Seddegh et al., 2017b) revealed that during charging the buoyancy forces induce the liquid PCM to move upwards; therefore, the upper section of the storage is affected. The concept of the method is to attribute the deviation between the two models (i.e. CCNC and PC) to natural convection:
γConv = γCCNC − γPC
γ=
γ U + 2γ L 3
(5)
Fig. 10 presents the flowchart followed for method development under SCD.
(4)
5.1.2. Verification of logic In this section, melting front separation for U and L (as shown in Fig. 9) is examined to check whether it would make reasonable sense in reality. To do so, some scenarios were selected as shown in Table 8. Then, the results obtained from the models were examined to verify how realistic such logic is. In Fig. 11, the right and left sides of each contour show the results for Case 16 (as an example) obtained from the
5.1.1. Logic description When SCD is initiated, heat transfer is dominated by conduction. However, as the PCM melts, convection gradually establishes, dominating the phase change mechanism (Joybari et al., 2017c). Note that under SCD, due to the simultaneous cooling from the outer tube, some
Fig. 4. Contours of liquid fraction for Cases 4–6 using the PC model. 793
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
Fig. 5. Contour of liquid fraction for Cases 4–6 using the CCNC model.
PC and CCNC models, respectively. There is an acceptable agreement between the reality and the simple method. Quantitative results are shown in Fig. 12. It should be noted that unlike the existing method (see Fig. 8), the new method is more robust especially at the early stages of SCD.
R=
(7)
Regression analysis (multiple variable) by least-squares method was applied to the data for liquid fraction of U. Fig. 14 shows the results with confidence level of 95% while the obtained constants for Eq. (6) are presented in Table 9. According to Fig. 14, the correlation can predict about 90% of data within a discrepancy of ± 15%.
5.2. Correlation development
5.2.2. Correlation for L Correlation development for L is more complicated than U. The reason is that some modifications are needed to come up with a similar correlation format presented in Eq. (6) for L. Fig. 15 shows an example (for Case 17) of why modifications are needed. According to Fig. 15, since the correlation for L (i.e. f2) initiates once the melting in U has completed (see the vertical line in green), the liquid fraction of the PC model would be equal to the last of the U correlation (i.e. f1). In other words, each case has its unique value for U ,f liquid fraction of the PC model (denoted by γPC ), achieving which U reaches complete melting. This point would be the initial point of the correlation for L. To be consistent with the correlation for U, the origin’s coordinate should be relocated to that unique point (see the red dot in the figure). Therefore, for all the cases the initial liquid fraction value would be zero:
The results obtained for each part of the storage are first presented, based on which correlation development is carried out.
5.2.1. Correlation for U The liquid fraction values of U are plotted versus the PC model in Fig. 13. Note that for some cases the maximum liquid fraction value is lower than one. This is due to a modification according to the behavior of heat transfer under SCD. As the propagation of melting front continues, the resistance against heat transfer from the HTF tube to the melting front increases, decreasing the rate of change for liquid fraction over time (due to reaching the steady state). Consequently, after a while, the steepness of the graph slope would increase. The moment for such behavior is affected by the PCM thermophysical properties, the storage geometrical properties as well as the considered boundary conditions. For the sake of method simplicity, carefully examining the data revealed that if the rate of change in liquid fraction of the PC model is beyond 10-5, the correlations can be developed. Although Fig. 13 shows the results for different geometry/PCM combinations, only three lines are distinguishable. Each of the lines represent a geometry shown in Table 2. This confirms only geometrical parameters affect the liquid fraction of U when plotted versus the PC model. Since the effect of PCM thermophysical properties are negligible, the correlation format can be:
γ U = a (R)b (γPC )c
ro ri
U ,f mod γPC = γPC − γPC
(8)
U ,f γ L, mod = γ L − γPC
(9)
For the modified L, Fig. 16 shows the results for different geometry/ PCM combinations, where the lines are hardly distinguishable due to reaching the steady state. Therefore: mod c γ L, mod = a (R)b (γPC )
(10)
The results of regression are shown in Fig. 17 and Table 10 tabulates the constant values of Eq. (10). According to Fig. 17, the correlation can predict all the data within a discrepancy of ± 3%. In Fig. 16, due to reaching the steady state under SCD, the data points do not exceed a certain modified PC model liquid fraction (i.e.
(6)
where the constants (i.e. a, b and c) should be obtained from regression. The storage radius ratio (R) is: 794
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
Fig. 8. Logic verification results of each scenario in Table 7 for Case 16 using the existing melting method.
Fig. 9. (a) The concept of two separate fronts for SCD in a TTHX and (b) the propagation of melting front before complete melting of U.
words, the presented simple front tracking model is a trade-off between simplicity and accuracy where the correlation for U (i.e. f1) can represent the melting front location almost alone, without that of L (i.e. f2).
Fig. 6. The existing procedure to develop the method of front tracking under melting for STHXs equipped with PCM (Joybari et al., 2017b).
Table 7 Logic verification scenarios for the existing melting method.
5.3. Verification
Scenario
Upper
Lower
Scenario 1 Scenario 2 Scenario 3
0.5 0.75 0.83
0.1412 0.1617 0.1712
So far three PCMs (shown in Table 3) where used for correlation development. In this section, a new PCM (i.e. RT35HC) is used to verify the correlations. The verification results for the three verification cases (shown by “Ver” in Table 4) are presented in Fig. 18 for U. It shows that the correlation can predict almost all results with a discrepancy of ± 15%. The verification results for L are shown in Fig. 19. Note that two cases (i.e. Ver 1 and Ver 3 in Table 4) achieved steady state earlier in Fig. 18. Therefore, they are missing points in Fig. 19. It can be concluded that the correlation can excellently predict the results with a discrepancy of ± 3%.
the x-axis) value of 0.05. Consequently, since the natural convection is decelerated, the prediction accuracy is greatly improved in Fig. 17. This is very important as it confirms that the sectionalizing of U and L as presented in this paper almost coincides with the steady state. In other
Fig. 7. Logic verification contours of each scenario in Table 7 for Case 16 using the existing melting method. 795
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
Fig. 12. Logic verification results of each scenario in Table 8 for Case 16 using the SCD method.
Fig. 13. Liquid fraction variation of U versus the PC model.
6. Conclusion Fig. 10. The procedure followed for method development under SCD.
This paper examined the applicability of an existing method for front tracking of PCMs under SCD. The performance of the existing method was found to improve towards the steady state condition, while SCD does not normally continue for such long periods. As a result, a new method for tracking the melting front was presented according to the natural convection behavior during SCD. The main findings/ achievements of this paper include:
Table 8 Logic verification scenarios for the SCD method. Scenario
U
L
Scenario 1 Scenario 2 Scenario 3
0.5 1 1
0.1285 0.1589 0.2523
• Separation of melting front should be carried out according to the natural convection to be more compatible with the SCD process. • Correlations were developed and verified for different portions of the storage (i.e. U and L) in a horizontal system. • The PCM thermophysical properties were found to be insignificant and only the storage geometrical properties were considered in the
Fig. 11. Logic verification contours of each scenario in Table 8 for Case 16 using the SCD method. 796
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
Fig. 14. Liquid fraction variation of U versus the correlation.
Fig. 17. Liquid fraction variation of modified L versus the correlation.
Table 9 Constants of Eq. (6) for the liquid fraction of U.
Table 10 Constants of Eq. (10) for the liquid fraction of modified L.
a
b
c
R2
Range
a
b
c
R2
Range
4.031
3.574
2.410
0.898
R = 1.19 – 2.31
1.580
1.846
0.979
0.9626
R = 1.19 – 2.31
Fig. 18. Verification of predictions for the correlation developed for U (i.e. f1). Fig. 15. Variation of liquid fractions for an exemplary case (i.e. Case 17).
Fig. 19. Verification of predictions for the correlation developed for modified L (i.e. f2).
Fig. 16. Liquid fraction variation of U versus the PC model (modified).
797
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
correlations (i.e. PC model liquid fraction and storage radius ratio).
energy storage system. Renew. Energy. Katsman, L., Dubovsky, V., Ziskind, G., Letan, R., 2007. Experimental Investigation of Solid-Liquid Phase Change in Cylindrical Geometry. In: ASME/JSME 2007 Thermal Engineering Heat Transfer Summer Conference collocated with the ASME 2007 InterPACK Conference, American Society of Mechanical Engineers, pp. 239–244. Klimeš, L., Mauder, T., Charvát, P., Štětina, J., 2018. Front tracking in modelling of latent heat thermal energy storage: assessment of accuracy and efficiency, benchmarking and GPU-based acceleration. Energy 155, 297–311. Li, C.-Y., Garimella, S.V., Simpson, J.E., 2003. Fixed-grid front-tracking algorithm for solidification problems, part I: method and validation. Numer. Heat Transf., Part B: Fundam. 43 (2), 117–141. Liu, C., Groulx, D., 2014. Experimental study of the phase change heat transfer inside a horizontal cylindrical latent heat energy storage system. Int. J. Therm. Sci. 82, 100–110. Liu, Z., Wang, Z., Ma, C., 2006. An experimental study on the heat transfer characteristics of a heat pipe heat exchanger with latent heat storage Part II: simultaneous charging/ discharging modes. Energy Convers. Manage. 47 (7), 967–991. Liu, Z., Yu, Z.J., Yang, T., Qin, D., Li, S., Zhang, G., Haghighat, F., Joybari, M.M., 2018. A review on macro-encapsulated phase change material for building envelope applications. Build. Environ. Long, J.-Y., Zhu, D.-S., 2008. Numerical and experimental study on heat pump water heater with PCM for thermal storage. Energy Build. 40 (4), 666–672. Mahdi, J.M., Nsofor, E.C., 2018. Solidification enhancement of PCM in a triplex-tube thermal energy storage system with nanoparticles and fins. Appl. Energy 211, 975–986. Mahdi, J.M., Lohrasbi, S., Ganji, D.D., Nsofor, E.C., 2019. Simultaneous energy storage and recovery in the triplex-tube heat exchanger with PCM, copper fins and Al2O3 nanoparticles. Energy Convers. Manage. 180, 949–961. Mat, S., Al-Abidi, A.A., Sopian, K., Sulaiman, M., Mohammad, A.T., 2013. Enhance heat transfer for PCM melting in triplex tube with internal–external fins. Energy Convers. Manage. 74, 223–236. Mirzaei, P.A., Haghighat, F., 2012. Modeling of phase change materials for applications in whole building simulation. Renew. Sustain. Energy Rev. 16 (7), 5355–5362. Murray, R.E., Groulx, D., 2014. Experimental study of the phase change and energy characteristics inside a cylindrical latent heat energy storage system: Part 2 simultaneous charging and discharging. Renew. Energy 63, 724–734. Najafian, A., Haghighat, F., Moreau, A., 2015. Integration of PCM in domestic hot water tanks: optimization for shifting peak demand. Energy Build. 106, 59–64. Omojaro, P., Breitkopf, C., 2014. Investigating and modeling of simultaneous charging and discharging of a PCM heat exchanger. Energy Procedia 48, 413–422. Omojaro, A.P., Breitkopf, C., 2017. Study on solid liquid interface heat transfer of PCM under simultaneous charging and discharging (SCD) in horizontal cylinder annulus. Heat Mass Transf. 1–18. Rathod, M.K., Jyotirmay, B., 2015. Development of correlation for melting time of phase change material in latent heat storage unit. Energy Procedia 75, 2125–2130. Riley, D., Smith, F., Poots, G., 1974. The inward solidification of spheres and circular cylinders. Int. J. Heat Mass Transf. 17 (12), 1507–1516. Seddegh, S., Wang, X., Joybari, M.M., Haghighat, F., 2017a. Investigation of the effect of geometric and operating parameters on thermal behavior of vertical shell-and-tube latent heat energy storage systems. Energy 137, 69–82. Seddegh, S., Wang, X., Henderson, A.D., 2015. Numerical investigation of heat transfer mechanism in a vertical shell and tube latent heat energy storage system. Appl. Therm. Eng. 87, 698–706. Seddegh, S., Wang, X., Henderson, A.D., 2016. A comparative study of thermal behaviour of a horizontal and vertical shell-and-tube energy storage using phase change materials. Appl. Therm. Eng. 93, 348–358. Seddegh, S., Joybari, M.M., Wang, X., Haghighat, F., 2017b. Experimental and numerical characterization of natural convection in a vertical shell-and-tube latent thermal energy storage system. Sustain. Cities Soc. 35, 13–24. Sharifi, N., Faghri, A., Bergman, T.L., Andraka, C.E., 2015. Simulation of heat pipe-assisted latent heat thermal energy storage with simultaneous charging and discharging. Int. J. Heat Mass Transf. 80, 170–179. Solomon, A., 1979. Melt time and heat flux for a simple PCM body. Sol. Energy 22 (3), 251–257. Sparrow, E., Broadbent, J., 1982. Inward melting in a vertical tube which allows free expansion of the phase-change medium. ASME J. Heat Transf. 104 (2), 309–315. Stritih, U., Tyagi, V.V., Stropnik, R., Paksoy, H., Haghighat, F., Mastani Joybari, M., 2018. Integration of passive PCM technologies for net-zero energy buildings. Sustain. Cities Soc. 41, 286–295. Tay, N., Bruno, F., Belusko, M., 2012c. Experimental validation of a CFD and an ε-NTU model for a large tube-in-tank PCM system. Int. J. Heat Mass Transf. 55 (21), 5931–5940. Tay, N., Belusko, M., Bruno, F., 2012a. An effectiveness-NTU technique for characterising tube-in-tank phase change thermal energy storage systems. Appl. Energy 91 (1), 309–319. Tay, N., Belusko, M., Bruno, F., 2012b. Experimental investigation of tubes in a phase change thermal energy storage system. Appl. Energy 90 (1), 288–297. Voller, V.R., Prakash, C., 1987. A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems. Int. J. Heat Mass Transf. 30 (8), 1709–1719. Yang, J., Yang, L., Xu, C., Du, X., 2016. Experimental study on enhancement of thermal energy storage with phase-change material. Appl. Energy 169, 164–176.
• The method can predict the location of the front under SCD outperforming the effective thermal conductivity method. Therefore, the developed method can replace the effective thermal conductivity method especially due to its simplicity. Overall, the users can obtain results within a discrepancy of 15%.
Acknowledgement The authors would like to express their gratitude to Concordia University for the support through Concordia Research Chair – Energy and Environment. References Abdulateef, A.M., Abdulateef, J., Mat, S., Sopian, K., Elhub, B., Mussa, M.A., 2018. Experimental and numerical study of solidifying phase-change material in a triplextube heat exchanger with longitudinal/triangular fins. Int. Commun. Heat Mass Transf. 90, 73–84. Al-abidi, A.A., Mat, S.B., Sopian, K., Sulaiman, M., Mohammed, A.T., 2013a. CFD applications for latent heat thermal energy storage: a review. Renew. Sustain. Energy Rev. 20, 353–363. Al-Abidi, A.A., Mat, S., Sopian, K., Sulaiman, M., Mohammad, A.T., 2013b. Experimental Investigation of Melting in Triplex Tube Thermal Energy Storage. Latest Trends in Renewable Energy and Environmental Informatics. Al-Abidi, A.A., Mat, S., Sopian, K., Sulaiman, M., Mohammad, A.T., 2013c. Experimental study of PCM melting in triplex tube thermal energy storage for liquid desiccant air conditioning system. Energy Build. 60, 270–279. Al-Abidi, A.A., Mat, S., Sopian, K., Sulaiman, M., Mohammad, A.T., 2013d. Internal and external fin heat transfer enhancement technique for latent heat thermal energy storage in triplex tube heat exchangers. Appl. Therm. Eng. 53 (1), 147–156. Al-Abidi, A.A., Mat, S., Sopian, K., Sulaiman, M., Mohammad, A.T., 2013e. Numerical study of PCM solidification in a triplex tube heat exchanger with internal and external fins. Int. J. Heat Mass Transf. 61, 684–695. Al-Abidi, A.A., Mat, S., Sopian, K., Sulaiman, M., Mohammad, A.T., 2014. Experimental study of melting and solidification of PCM in a triplex tube heat exchanger with fins. Energy Build. 68, 33–41. ANSYS Fluent Theory Guide, ANSYS, Inc., Canonsburg, PA, USA, 2013. Bruno, M.B.F., 2008. Design methodology of PCM thermal storage systems with parallel plates. EUROSUN 2008. El-Sawi, A., Haghighat, F., Akbari, H., 2014. Assessing long-term performance of centralized thermal energy storage system. Appl. Therm. Eng. 62 (2), 313–321. Fan, L.-W., Zhu, Z.-Q., Zeng, Y., Lu, Q., Yu, Z.-T., 2014. Heat transfer during melting of graphene-based composite phase change materials heated from below. Int. J. Heat Mass Transf. 79, 94–104. Rubitherm Technologies GmbH, in: Organic PCM, http://www.rubitherm.eu/en/index. php/productcategory/organische-pcm-rt/, Berlin, Germany, 2017. The Copper Tube Handbook, Copper Development Association Inc., New York, NY, USA, 2011. Hasan, A., 1994a. Thermal energy storage system with stearic acid as phase change material. Energy Convers. Manage. 35 (10), 843–856. Hasan, A., 1994b. Phase change material energy storage system employing palmitic acid. Sol. Energy 52 (2), 143–154. Ismail, K.A., Maria das Graças, E., 2003. Melting of PCM around a horizontal cylinder with constant surface temperature. Int. J. Therm. Sci. 42 (12), 1145–1152. Jian-you, L., 2008. Numerical and experimental investigation for heat transfer in triplex concentric tube with phase change material for thermal energy storage. Sol. Energy 82 (11), 977–985. Jones, B.J., Sun, D., Krishnan, S., Garimella, S.V., 2006. Experimental and numerical study of melting in a cylinder. Int. J. Heat Mass Transf. 49 (15), 2724–2738. Joybari, M.M., Haghighat, F., Moffat, J., Sra, P., 2015. Heat and cold storage using phase change materials in domestic refrigeration systems: the state-of-the-art review. Energy Build. 106, 111–124. Joybari, M.M., Haghighat, F., 2016. Natural convection modelling of melting in phase change materials: Development of a novel zonal effectiveness-NTU. 9th International Conference on Indoor Air Quality Ventilation and Energy Conservation In Buildings (IAQVEC 2016), Incheon Songdo, South Korea. Joybari, M.M., Haghighat, F., Seddegh, S., Al-Abidi, A.A., 2017a. Heat transfer enhancement of phase change materials by fins under simultaneous charging and discharging. Energy Convers. Manage. 152, 136–156. Joybari, M.M., Haghighat, F., Seddegh, S., 2017b. Natural convection characterization during melting of phase change materials: development of a simplified front tracking method. Sol. Energy 158, 711–720. Joybari, M.M., Haghighat, F., Seddegh, S., 2017c. Numerical investigation of a triplex tube heat exchanger with phase change material: simultaneous charging and discharging. Energy Build. 139, 426–438. Joybari, M.M., Seddegh, S., Wang, X., Haghighat, F., 2019. Experimental investigation of multiple tube heat transfer enhancement in a vertical cylindrical latent heat thermal
798
Contents lists available at ScienceDirect
Solar Energy journal homepage: www.elsevier.com/locate/solener
Simultaneous charging and discharging of phase change materials: Development of correlation for liquid fraction
T
⁎
Mahmood Mastani Joybaria, Fariborz Haghighata, , Saeid Seddeghb, Yanping Yuanc a
Energy and Environment Group, Department of Building, Civil and Environmental Engineering, Concordia University, Montreal H3G 1M8, Canada School of Engineering, College of Science and Engineering, University of Tasmania, Hobart, TAS 7001, Australia c School of Mechanical Engineering, Southwest Jiaotong University, Chengdu 610031, China b
A R T I C LE I N FO
A B S T R A C T
Keywords: Front tracking Phase change material Natural convection Triplex tube heat exchanger Simultaneous charging and discharging
In this study, a method to track the melting front of phase change materials (PCMs) is developed under simultaneous charging and discharging (SCD) inside horizontal triplex tube heat exchangers. For SCD, the heat transfer mode of internal heating and external cooling is investigated. According to the heat transfer mechanism within the PCMs, two models of pure conduction (PC) and combined conduction and natural convection (CCNC) are developed and solved using ANSYS Fluent (version 17.2). Generally, during charging (i.e. melting), the buoyancy forces induce the liquid PCM to move upwards; therefore, the upper section of the storage is affected. However, this motion is limited during SCD due to the simultaneous cooling from the outer tube. Therefore, a front tracking method is developed assuming that natural convection only contributes melting to the upper section of the storage unit. On the other hand, it is assumed that the melting front for the lower section behaves similar to the PC model. In this study, three PCMs are used in three geometries to develop the correlations using two dimensionless numbers; i.e. liquid fraction of the PC model and the storage radius ratio. Another PCM is then used to verify the method. The correlations can provide results in the range of ± 15% discrepancy.
1. Introduction Latent heat storage using phase change materials (PCMs) benefits from high storage capacity and can adjust the time mismatch between energy supply and demand (Mirzaei and Haghighat, 2012). Therefore, it found several applications, including building envelopes (Liu et al., 2018), ventilation systems (El-Sawi et al., 2014), refrigeration systems (Joybari et al., 2015), net zero energy buildings (Stritih et al., 2018), hot water tanks (Najafian et al., 2015), etc. Prior to designing any latent heat storage system using PCMs, the heat transfer mechanisms should be understood. In the literature, heat transfer within PCMs can be modeled using (i) the pure conduction (PC) model or (ii) combined conduction and natural convection (CCNC) model. The PC model was used in earlier studies to investigate the heat transfer within PCMs (Seddegh et al., 2015). Since gravity has no effect on conduction, the melting front would be circular (or cylindrical in 3D) (Tay et al., 2012a; Joybari and Haghighat, 2016). In reality, however, as the PCM melts, buoyancy forces due to the change in PCM density create an upward motion (Seddegh et al., 2016). This
phenomenon is known as natural convection, which has been observed in several experimental studies (Al-Abidi et al., 2013b; Joybari et al., 2019). The two main design tool categories for heat exchangers are (i) simplified methods, and (ii) complex methods including CFD simulations. The first category includes the well-known approaches such as effectiveness-number of transfer units (ε-NTU), normally oversimplifying the heat transfer process within PCMs. For instance, neglecting natural convection results in effectiveness underestimation during the phase change process (Tay et al., 2012a,b,c; Bruno, 2008). In the second category, CFD simulations are more complicated and normally require long computational time. Nevertheless, some tools can be developed, whose simplicity and accuracy lie between the two mentioned categories. To enhance the accuracy of the ε-NTU method, a method called effective thermal conductivity was introduced which accounts for the natural convection during the melting process. In this method, the thermal conductivity values are artificially increased to mimic the experimental results. Comparison of ε-NTU method with experimental
Abbreviations: CFD, computational fluid dynamics; CCNC, combined conduction and natural convection; HTF, heat transfer fluid; NTU, number of transfer units; PC, pure conduction; PCM, phase change material; STHX, shell and tune heat exchanger; TTHX, triplex tube heat exchanger ⁎ Corresponding author. E-mail address: [email protected] (F. Haghighat). https://doi.org/10.1016/j.solener.2019.06.051 Received 6 March 2019; Received in revised form 28 May 2019; Accepted 20 June 2019 0038-092X/ © 2019 International Solar Energy Society. Published by Elsevier Ltd. All rights reserved.
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
density (kg m−3)
Nomenclature
ρ
a b c C∗ Fo h∗ P* Pr r R Ra Ste V
Subscripts
constant constant constant dimensionless Darcy coefficient (-) Fourier number (-) dimensionless enthalpy (-) dimensionless pressure (-) Prandtl number (-) radius (m) storage radius ratio (-) Rayleigh number (-) Stefan number (-) dimensionless velocity (-)
Superscripts final L section modified U section dimensionless variable
f L mod U ∗
Greek symbols
γ
combined conduction and natural convection model inner outer pure conduction model reference
CCNC i o PC ref
liquid fraction (-)
In this method, both phases (i.e. solid and liquid) are modeled using the same enthalpy-based energy equation (Voller and Prakash, 1987). Therefore, to track the melting front using the enthalpy method, a twostep process should be followed: (i) calculation of the enthalpy values, and (ii) determination of the melting front location according to the temperature values (Li et al., 2003). Nevertheless, despite the progress in computational power, the numerical analysis of PCMs is still complicated. It also requires code development expertise or accessibility to a commercial software for simulation. Therefore, to characterize the PCM phase change process, several correlations have been reported in the literature (see Table 1). As the table indicates, the older correlations were developed for Fourier number (Fo) to obtain the phase change duration. However, more recent studies focused on liquid fraction (γ). Note that natural convection was not considered in most of the studies. Overall, although correlations have been developed for PC and CCNC models, to the best of the authors’ knowledge, no correlation has been developed linking these two models. Another inference from Table 1 is that PCMs have normally been investigated under one-time (or consecutive processes of) charging and/or discharging. In other words, an initially solidified (or melted) PCM was considered to undergo melting (or solidification) process. This ideal assumption is not in line with real implementations where the PCM might undergo simultaneous charging and discharging (SCD). In such cases, a PCM is being melted from one side, while it is
data revealed that the method would have acceptable results only if it is enhanced by the effective thermal conductivity method (Tay et al., 2012c). Nevertheless, there are some drawbacks for the effective thermal conductivity method: (i) it requires experimental data; therefore, an experimental setup should be configured and manufactured, (ii) it is time-consuming as the best value is obtained by trial and error while comparing to experimental data, and (iii) despite all these complexities, since the effective thermal conductivity method is based on the conduction concept, no information can be obtained regarding the melting front. Knowing the location of melting front can not only indicate the liquid fraction value, but also capture partial melting. However, obtaining liquid fraction in experimental studies is very complicated. During earlier experimental studies the phase change process was interrupted after some time in order to extract the remaining solid part (Sparrow and Broadbent, 1982). In recent years, transparent tubes were utilized which can provide direct visual observation (Liu and Groulx, 2014). Moreover, digital high resolution photography (Yang et al., 2016) as well as image processing (Jones et al., 2006) has greatly facilitated the liquid fraction calculation process. Nevertheless, it is preferable to conduct numerical analyses in order to obtain liquid fraction. Melting front tracking used to be a great challenge in early numerical studies since a moving-boundary problem (due to the change of melting front location) used to be solved (Klimeš et al., 2018). However, such a complexity can be avoided by using the enthalpy method.
Table 1 Hierarchical summary of the existing correlations characterizing the phase change process. Refs.
Type of study
Storage type(s)
Orientation
(Riley et al., 1974)
Analytical
CylindricalSpherical
Fo = a +
b Ste
(Solomon, 1979)
Analytical
Fo = a +
b Ste
(Hasan, 1994b)
Experimental
Slab CylindricalSpherical Cylindrical
Fo = a +
b Ste
(Hasan, 1994a)
Experimental
Cylindrical
Fo = a +
b Ste
(Ismail and Maria das Graças, 2003)
Numerical
Cylindrical
Horizontal Vertical Horizontal Vertical Horizontal
(Katsman et al., 2007)
ExperimentalNumerical
STHX
Vertical
γ = (FoSte )aRaDb exp
γ = aFoSte bRac Fo = a + b Re + cSte
(Fan et al., 2014)
Experimental
Cylindrical
Vertical
(Rathod and Jyotirmay, 2015)
Experimental
STHX
Vertical
789
Correlation
Comment Developed for discharging (no significant natural convection occurred). Developed for both charging and discharging.
Y = a (X )2 + b (X ) + c
( ) c RaH
Developed for both charging and discharging. Natural convection was neglected despite reporting its effect. Developed for both charging and discharging. Natural convection was neglected despite reporting its effect. Developed for charging. The correlations include complete melting duration and melted volume Stefan (Ste) or Rayleigh (Ra) number is represented by X. Developed for charging. Developed for charging. Developed for charging.
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
• No correlation has been developed for PCMs linking the CCNC and
simultaneously being solidified from the other. This can happen when a PCM is unpredictably and/or intermittently charged and/or discharged. As an example, when PCM is used for solar domestic hot water, occupants might demand water consumption (i.e. discharging) during sunshine (i.e. charging) hours. According to the literature, despite several studies which investigated the charging, discharging or consecutive charging and discharging process, the behavior of PCMs when undergoing SCD has not been thoroughly studied. For instance, Liu et al. (2006) experimentally investigated the SCD process in a heat pipe heat exchanger. A paraffin PCM was used where its high thermal resistance was found to inhibit efficient system performance. However, that was not the case for a metal high-temperature PCM in another study (Sharifi et al., 2015). Murray and Groulx (2014) conducted experimental investigations for SCD of a PCM for domestic hot water applications. Charging was continuously imposed, while the discharging was intermittent. Again, the solid part of the PCM controlled the heat exchange process due to its high thermal resistance. Omojaro and Breitkopf (2014) numerically and experimentally investigated SCD in a horizontal shell and tube heat exchanger (STHX) where the PCM was located in the shell. For charging and discharging, a heater and a fan were used, respectively. For numerical investigation, a simple symmetric one-dimensional model was developed for temperature as well as liquid fraction. Recently, the same authors used the existing experimental setup to investigate the system performance under different convection mechanisms from the shell (i.e. natural convection and forced convection) (Omojaro and Breitkopf, 2017). It was reported that forced convection resulted in lower melting rate. Besides, lower air velocities caused shorter transition time within the PCM from conduction to convection. Thus, the low thermal conductivity of the PCM deteriorated the system performance. Certain configurations have the capability to undergo SCD. Triplex tube heat exchangers (TTHXs) consist of three concentric tubes. For heat storage, the PCM fills the middle space, while the HTFs flow through the inner and outer tubes. TTHXs also possess some characteristics, suitable for solar domestic hot water applications: (i) possibility to impose SCD, (ii) possessing separate tubes for the HTF and water to flow which can satisfy sanitary requirements, and (iii) enhanced natural convection, improving heat transfer within the solid PCM. Several studies investigated the phase change process in PCMequipped TTHXs under one-time solidification (Al-Abidi et al., 2013e; Abdulateef et al., 2018; Mahdi and Nsofor, 2018), one-time melting (AlAbidi et al., 2013b, 2013c, 2013d; Mat et al., 2013), or consecutive melting and solidification (Al-Abidi et al., 2014; Jian-you, 2008; Long and Zhu, 2008). In addition, SCD in TTHXs has recently gained attention. For instance, Mahdi et al. (2019) recently optimized the fin configuration in TTHXs under SCD and reported that it outperformed addition of nanoparticles (i.e. a PCM thermal conductivity enhancement technique). In our earlier study, a horizontal TTHX was numerically investigated under SCD (Joybari et al., 2017c). For internal heating/ external cooling, it was found that the PCM can be partially melted in the upper half of the system. Later, fin configurations were recommended under SCD as well as non-SCD applications for TTHXs in accordance with the behavior of natural convection (Joybari et al., 2017a). Reviewing the available literature demonstrates the following limitations:
• • • •
PC models under SCD.
Therefore, in this study, a simple method of front tracking is developed for SCD, considering the natural convection behavior. In our earlier study (Joybari et al., 2017b), a method was developed for melting process only. Therefore, this study extends its application to SCD providing important information regarding the melting front location and partial melting. 2. Materials and methods The considered storage geometrical properties as well as PCM thermophysical properties are described in this section, which are the same as our earlier study for melting (Joybari et al., 2017b). A summary of the modeling approach is also presented. 2.1. Setup characteristics The US standard for sizing of copper tubes was used which recommends tubes of “Type M” for water heating applications (Handbook et al., 2011). A horizontal TTHX was investigated in which the PCM filled the middle tube (see Fig. 1). Three geometries (see Table 2) were used for data generation (to be used for correlation development). Note that the inner tube radius (ri) was altered, while the other tube size (ro) remained the same. 2.2. PCM thermophysical properties In this study, four commercially available paraffin PCMs (from Rubitherm Technologies GmbH) were investigated (see Table 3). Note that during data generation process the highlighted PCM (i.e. RT35HC) was excluded so that it can be used later for correlation verification. 2.3. Governing equations This section discusses the assumptions as well as the modeling procedure adopted for this study. 2.3.1. Assumptions The following were considered for the phase change process of the PCMs:
• Assumed
laminar, incompressible Newtonian fluid flow for the
ri
Most earlier studies (especially older ones) neglected the natural convection during charging and/or discharging process, CFD simulations using the CCNC model are not only time consuming but also require high computational power, The phase change process under SCD has not been thoroughly studied in the literature and experimental or numerical investigations of PCM under SCD are limited, especially in TTHXs, There is no correlation in the literature for front tracking under SCD,
ro ș=0
Fig. 1. Schematic representation of the TTHXs. 790
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
Table 2 Geometrical properties of the middle space in the TTHXs (Joybari et al., 2017b; Handbook et al., 2011). Geometry
Geometry 1 Geometry 2 Geometry 3
Radius (mm)
Storage radius ratio (R)
Inner (ri)
Outer (ro)
26.98 39.69 52.38
62.32 62.32 62.32
2.31 1.57 1.19
melted PCM,
• Neglected viscous dissipation due to the absence of large velocity gradients, • Assumed no-slip and Dirichlet boundary conditions for the fluid flow and heat transfer, respectively, • Assumed linear variation of liquid fraction with temperature, • Assumed the thermophysical properties to be constant for each
Fig. 2. Model validation under the melting process.
solidification processes. The simulation results are compared with experimental data obtained from two earlier studies.
phase except density, being evaluated using Boussinesq approximation (considering natural convection).
2.3.4.1. Melting process. The experimental data of Mat et al. (2013) was used to validate the model under melting. The paraffin-based RT82 (Rubitherm Technologies GmbH) was used in that study as the PCM in a horizontal finned TTHX. To melt the PCM, hot water was used, flowing inside both inner and outer tubes. Experimental data (i.e. average temperature) is compared with numerical results in Fig. 2. Despite slight overestimation of the model, an acceptable agreement exists with the experimental data. Therefore, the model is considered validated under melting.
2.3.2. Modeling Our earlier studies (Joybari et al., 2017a,c) include the detailed modeling approach. For the sake of conciseness and generalization, this section presents the equations using dimensionless parameters: Continuity Equation:
→ ∂ρ + ∇ ·(ρV ) = 0 ∂Fo
(1)
Momentum Balance:
→ → →→ → → ∂V Ra ·Pr ∗ ∗ + ∇ ·( V V ) = −∇P ∗ + Pr ∇ ·(∇ V ) + (h − href ) − C ∗V ∂Fo Ste
2.3.4.2. Solidification process. Experimental data of Seddegh et al. (2017a) was used to validate the model under solidification. The paraffin-based RT60 (Rubitherm technologies GmbH) was used as the PCM in a vertical STHX. Water was used as the HTF in the tube to solidify the PCM. The temperature of an exact point (i.e. T8 of Cylinder D in Seddegh et al. (2017a) was considered for validation. Fig. 3 shows the results where the model followed the trend of experimental data. Note that the deviations are because any thermocouple records the temperature of its surrounding chunk of the PCM. In other words, they are not capable of measuring the point-specific temperature. Despite this, the numerical results reasonably agree with those obtained experimentally. Therefore, the model is considered validated under solidification as well.
(2)
Energy Balance: ∗ → ∂γ ⎡ ∂h + ∇ ·(→ + ∇ ·( V γ ) ⎤ = ∇ ·(∇h∗) V h∗) ⎤ + ⎡ ⎣ ∂Fo ⎦ ⎣ ∂Fo ⎦
(3)
2.3.3. Numerical approach Several commercial software exists for numerical investigation of the phase change process in PCMs. Since ANSYS Fluent (based on the enthalpy method) is most commonly used (Al-abidi et al., 2013a), it was adopted for this study. First, the 2D geometry was developed in DesignModeler. Thereafter, ANSYS Meshing was used to generate triangular meshes. Then, simulation setup was carried out in ANSYS Fluent. The pressure-based solver was utilized as it is compatible with solving phase change processes (ANSYS Fluent Theory Guide, 2013). For pressure-velocity coupling, the SIMPLE method was used. QUICK and PRESTO schemes were adopted to spatially discretize the momentum and pressure fields, respectively.
2.3.5. Analysis procedure In this study, the simulations were performed for each PCM (details in Table 3) in various geometries (details in Table 2). Consequently, Table 4 tabulates all resulting combinations for Stefan numbers as well as storage radius ratios. Note that in Table 4 three cases are designated as “Ver”. These indicate cases which are later used for correlation verification. Overall, 42 simulations were conducted as both the PC and CCNC models were simulated.
2.3.4. Model validation In this section, the model is validated separately for melting and
Table 3 The PCMs and their thermophysical properties (Joybari et al., 2017b; GmbH, 2017). Property
Dimension
RT31
RT35
RT35HC
RT44HC
Melting onset temperature Solidification onset temperature Solid density Liquid density Specific heat capacity Latent heat Thermal conductivity Thermal expansion coefficient Dynamic viscosity
K (°C) K (°C) kg/m3 kg/m3 J/kg K J/kg W/m K 1/K kg/m s
300.15 (27) 306.15 (33) 880 760 2000 170,000 0.2 0.00076 0.002508
302.15 (29) 309.15 (36) 860 770 2000 160,000 0.2 0.00076 0.002500
307.15 (34) 309.15 (36) 880 770 2000 240,000 0.2 0.00076 0.002700
314.15 (41) 317.15 (44) 800 700 2000 250,000 0.2 0.00076 0.003300
791
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
cooling from the outer tube, maintaining a solid PCM layer adjacent to the outer tube wall. 3.3. Comparison According to the results, compared to the CCNC model, the PC model benefits from:
• Easier analysis setup as well as shorter running time, • Easier meeting of convergence criteria for its energy equation (since continuity as well as momentum equations are not present), • Possibility of analytical solution (under certain circumstances), drastically reducing the computational time.
According to these, the present study is focused on developing a simple method, which can generate the CCNC model results using PC model data. The method can at the same time approximately locate the melting front, outperforming the effective thermal conductivity method.
Fig. 3. Model validation under the solidification process. Table 4 The considered cases indicating their Stefan number and storage radius ratio. Case
R
Ste
Case
R
Ste
Case
R
Ste
Case 1 Case 2 Case 3 Case 4 Case 5 Case 6 Ver 1
2.31 1.57 1.19 2.31 1.57 1.19 2.31
0.16 0.16 0.16 0.24 0.24 0.24 0.33
Case 7 Case 8 Case 9 Case 10 Case 11 Case 12 Ver 2
2.31 1.57 1.19 2.31 1.57 1.19 1.57
0.25 0.25 0.25 0.32 0.32 0.32 0.17
Case 13 Case 14 Case 15 Case 16 Case 17 Case 18 Ver 3
2.31 1.57 1.19 2.31 1.57 1.19 1.19
0.47 0.47 0.47 0.50 0.50 0.50 0.33
4. Performance analysis of the existing method In our earlier study (Joybari et al., 2017b), a method of front tracking was developed for melting in STHXs. The same geometry configurations as this study were used to develop the correlations. Fig. 6 shows the procedure used in that study for the method development. In this section, the possibility for application of the existing method under SCD is verified. According to Fig. 6, the existing method separated the melting front into two parts for each half of the storage based on the mechanism of natural convection. Since natural convection affects the upper half of horizontal storages, its effect was first assumed to be solely on the upper half. Consequently, the behavior of the lower half was considered identical to the PC model. However, after the upper half was melted completely, the existing method deviated the lower half by attributing the remainder of natural convection. To verify the performance of the existing method under SCD, three scenarios have been considered as shown in Table 7. The first two scenarios consider liquid fractions of 0.5 and 0.75 for the upper half, while Scenario 3 is for the steady state condition under SCD. In Fig. 7 the right and left sides of each contour show the results for the PC and CCNC models, respectively. According to the contours, since the upper half of the system never reaches complete melting under SCD (as shown in Fig. 5), the existing method (Joybari et al., 2017b) cannot accurately capture the mechanism of natural convection. Moreover, Fig. 8 presents the quantitative comparison between the existing method and the observed one. As the figure shows, the performance of the existing method improves towards the steady state condition (i.e. Scenario 3) mainly due to the deceleration of natural convection. Note that a thermal energy storage system should not always run under SCD. The reason is that under such conditions there would be no benefit for the PCM. In other words, the HTFs can have direct heat transfer. Consequently, an ideal method should be capable of providing more accurate results specially at the early stages of SCD (as the steady state condition may never occur). Therefore, the existing method may be applied under SCD only if steady state condition is to be reached; otherwise, it should not
Before conducting the analysis, the simulation independency in terms of grid size and timestep size were checked. The results for Case 4 are presented in this section as an example. To check grid size independency three resolutions (i.e. 2574, 4091 and 7930 elements) were considered. In addition, three timestep sizes (i.e. 0.05, 0.1 and 0.2 s) were also investigated. Tables 5 and 6 tabulate the liquid fraction results for these criteria under SCD for 80 min. According to the results, 0.1 s timestep with 4091 elements was found to be sufficient for the simulation. 3. Results and discussion Prior to method description, its underlying concept is elaborated by explaining and comparing the mechanism of PC and CCNC models. 3.1. Pure conduction model As its name implies, the PC model only considers conduction and neglects natural convection during melting. Although the PC model is not realistic, its simplicity and fast solution (under certain circumstances it could be analytically solved) make it favorable. Fig. 4 illustrates the melting front propagation for RT31 under SCD using the PC model. The liquid fraction contours show that the melting front is circular as it monotonically propagates in the radial direction. 3.2. Combined conduction and natural convection model According to experimental data (Seddegh et al., 2017b) when a PCM melts, buoyancy forces are formed due to the density change, creating natural convection. Fig. 5 illustrates how the melting front propagates for RT31 under SCD using the CCNC model. Comparing Fig. 5 (for the CCNC model) with that of the PC model (i.e. Fig. 4) it can be understood that the upward liquid PCM motion affects the upper section of the storage, enhancing heat transfer and subsequently shortening the phase change process. Another very important observation is that despite the upward motion, for none of the cases the upper half achieved complete melting (i.e. a liquid fraction of one). This is due to the simultaneous
Table 5 Verification of grid size independency (with 0.1 s timestep). Number of elements
2574 4091 7930
792
Simulation time (s) 300
600
1200
2400
4800
0.0238 0.0278 0.0312
0.0442 0.0540 0.0589
0.0918 0.1098 0.1153
0.1927 0.2283 0.2202
0.3632 0.3777 0.3354
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
part of the PCM (adjacent to the outer tube) will remain solid; thus, complete melting cannot be achieved. A closer look at Fig. 5 reveals that within the upper half of the system, the effect of natural convection on the top part is stronger than the bottom part. To further clarify this, an exemplary contour is shown in Fig. 9a. Despite simultaneous cooling, the upward melted PCM motion could melt those parts of the PCM denoted by U (i.e. the upper one third of the storage). The melted region in U is greater than the rest (denoted by L). This type of sectionalizing has been defined to meet method simplicity requirements upon scrutinized examination of liquid fraction contours under SCD. Consequently, natural convection is first attributed to propagate the location of melting front at U (see Fig. 9b) until its liquid fraction reaches a value of one. In other words, there would be two separate melting fronts where the front radius of U is larger than that of L. The developed correlation for U (i.e. f1) would map the results of the PC model to those of U. Once the PCM in U reaches complete melting (i.e. a liquid fraction value of one), the remainder of natural convection is designated to L. In other words, L starts to deviate from the PC model once U is completely melted. The developed correlation for L (i.e. f2) would map the results of the PC model to those of L. The phase change process would later stabilize once the steady state is reached under SCD. For the two melting fronts shown in Fig. 9, the overall liquid fraction is:
Table 6 Verification of timestep size independency (with 4091 elements). Timestep (s)
0.05 0.1 0.2
Simulation time (s) 300
600
1200
2400
4800
0.0276 0.0278 0.0279
0.0537 0.0540 0.0542
0.1083 0.1098 0.1101
0.2228 0.2283 0.2302
0.3667 0.3777 0.3836
be used, as it is not compatible with natural convection under SCD. 5. Front tracking method development for simultaneous charging and discharging Due to the improper performance of the existing method, a new method of front tracking together with its liquid fraction correlations are developed in this section specifically for SCD. Thereafter, a new PCM is used for verification of the correlations. 5.1. Description of the method When natural convection is absent (i.e. the PC model), the shape of melting front would be cylindrical (i.e. circular in 2D). Consequently, any fraction of the storage would have the same liquid fraction as the PC model. However, experimental data (Seddegh et al., 2017b) revealed that during charging the buoyancy forces induce the liquid PCM to move upwards; therefore, the upper section of the storage is affected. The concept of the method is to attribute the deviation between the two models (i.e. CCNC and PC) to natural convection:
γConv = γCCNC − γPC
γ=
γ U + 2γ L 3
(5)
Fig. 10 presents the flowchart followed for method development under SCD.
(4)
5.1.2. Verification of logic In this section, melting front separation for U and L (as shown in Fig. 9) is examined to check whether it would make reasonable sense in reality. To do so, some scenarios were selected as shown in Table 8. Then, the results obtained from the models were examined to verify how realistic such logic is. In Fig. 11, the right and left sides of each contour show the results for Case 16 (as an example) obtained from the
5.1.1. Logic description When SCD is initiated, heat transfer is dominated by conduction. However, as the PCM melts, convection gradually establishes, dominating the phase change mechanism (Joybari et al., 2017c). Note that under SCD, due to the simultaneous cooling from the outer tube, some
Fig. 4. Contours of liquid fraction for Cases 4–6 using the PC model. 793
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
Fig. 5. Contour of liquid fraction for Cases 4–6 using the CCNC model.
PC and CCNC models, respectively. There is an acceptable agreement between the reality and the simple method. Quantitative results are shown in Fig. 12. It should be noted that unlike the existing method (see Fig. 8), the new method is more robust especially at the early stages of SCD.
R=
(7)
Regression analysis (multiple variable) by least-squares method was applied to the data for liquid fraction of U. Fig. 14 shows the results with confidence level of 95% while the obtained constants for Eq. (6) are presented in Table 9. According to Fig. 14, the correlation can predict about 90% of data within a discrepancy of ± 15%.
5.2. Correlation development
5.2.2. Correlation for L Correlation development for L is more complicated than U. The reason is that some modifications are needed to come up with a similar correlation format presented in Eq. (6) for L. Fig. 15 shows an example (for Case 17) of why modifications are needed. According to Fig. 15, since the correlation for L (i.e. f2) initiates once the melting in U has completed (see the vertical line in green), the liquid fraction of the PC model would be equal to the last of the U correlation (i.e. f1). In other words, each case has its unique value for U ,f liquid fraction of the PC model (denoted by γPC ), achieving which U reaches complete melting. This point would be the initial point of the correlation for L. To be consistent with the correlation for U, the origin’s coordinate should be relocated to that unique point (see the red dot in the figure). Therefore, for all the cases the initial liquid fraction value would be zero:
The results obtained for each part of the storage are first presented, based on which correlation development is carried out.
5.2.1. Correlation for U The liquid fraction values of U are plotted versus the PC model in Fig. 13. Note that for some cases the maximum liquid fraction value is lower than one. This is due to a modification according to the behavior of heat transfer under SCD. As the propagation of melting front continues, the resistance against heat transfer from the HTF tube to the melting front increases, decreasing the rate of change for liquid fraction over time (due to reaching the steady state). Consequently, after a while, the steepness of the graph slope would increase. The moment for such behavior is affected by the PCM thermophysical properties, the storage geometrical properties as well as the considered boundary conditions. For the sake of method simplicity, carefully examining the data revealed that if the rate of change in liquid fraction of the PC model is beyond 10-5, the correlations can be developed. Although Fig. 13 shows the results for different geometry/PCM combinations, only three lines are distinguishable. Each of the lines represent a geometry shown in Table 2. This confirms only geometrical parameters affect the liquid fraction of U when plotted versus the PC model. Since the effect of PCM thermophysical properties are negligible, the correlation format can be:
γ U = a (R)b (γPC )c
ro ri
U ,f mod γPC = γPC − γPC
(8)
U ,f γ L, mod = γ L − γPC
(9)
For the modified L, Fig. 16 shows the results for different geometry/ PCM combinations, where the lines are hardly distinguishable due to reaching the steady state. Therefore: mod c γ L, mod = a (R)b (γPC )
(10)
The results of regression are shown in Fig. 17 and Table 10 tabulates the constant values of Eq. (10). According to Fig. 17, the correlation can predict all the data within a discrepancy of ± 3%. In Fig. 16, due to reaching the steady state under SCD, the data points do not exceed a certain modified PC model liquid fraction (i.e.
(6)
where the constants (i.e. a, b and c) should be obtained from regression. The storage radius ratio (R) is: 794
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
Fig. 8. Logic verification results of each scenario in Table 7 for Case 16 using the existing melting method.
Fig. 9. (a) The concept of two separate fronts for SCD in a TTHX and (b) the propagation of melting front before complete melting of U.
words, the presented simple front tracking model is a trade-off between simplicity and accuracy where the correlation for U (i.e. f1) can represent the melting front location almost alone, without that of L (i.e. f2).
Fig. 6. The existing procedure to develop the method of front tracking under melting for STHXs equipped with PCM (Joybari et al., 2017b).
Table 7 Logic verification scenarios for the existing melting method.
5.3. Verification
Scenario
Upper
Lower
Scenario 1 Scenario 2 Scenario 3
0.5 0.75 0.83
0.1412 0.1617 0.1712
So far three PCMs (shown in Table 3) where used for correlation development. In this section, a new PCM (i.e. RT35HC) is used to verify the correlations. The verification results for the three verification cases (shown by “Ver” in Table 4) are presented in Fig. 18 for U. It shows that the correlation can predict almost all results with a discrepancy of ± 15%. The verification results for L are shown in Fig. 19. Note that two cases (i.e. Ver 1 and Ver 3 in Table 4) achieved steady state earlier in Fig. 18. Therefore, they are missing points in Fig. 19. It can be concluded that the correlation can excellently predict the results with a discrepancy of ± 3%.
the x-axis) value of 0.05. Consequently, since the natural convection is decelerated, the prediction accuracy is greatly improved in Fig. 17. This is very important as it confirms that the sectionalizing of U and L as presented in this paper almost coincides with the steady state. In other
Fig. 7. Logic verification contours of each scenario in Table 7 for Case 16 using the existing melting method. 795
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
Fig. 12. Logic verification results of each scenario in Table 8 for Case 16 using the SCD method.
Fig. 13. Liquid fraction variation of U versus the PC model.
6. Conclusion Fig. 10. The procedure followed for method development under SCD.
This paper examined the applicability of an existing method for front tracking of PCMs under SCD. The performance of the existing method was found to improve towards the steady state condition, while SCD does not normally continue for such long periods. As a result, a new method for tracking the melting front was presented according to the natural convection behavior during SCD. The main findings/ achievements of this paper include:
Table 8 Logic verification scenarios for the SCD method. Scenario
U
L
Scenario 1 Scenario 2 Scenario 3
0.5 1 1
0.1285 0.1589 0.2523
• Separation of melting front should be carried out according to the natural convection to be more compatible with the SCD process. • Correlations were developed and verified for different portions of the storage (i.e. U and L) in a horizontal system. • The PCM thermophysical properties were found to be insignificant and only the storage geometrical properties were considered in the
Fig. 11. Logic verification contours of each scenario in Table 8 for Case 16 using the SCD method. 796
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
Fig. 14. Liquid fraction variation of U versus the correlation.
Fig. 17. Liquid fraction variation of modified L versus the correlation.
Table 9 Constants of Eq. (6) for the liquid fraction of U.
Table 10 Constants of Eq. (10) for the liquid fraction of modified L.
a
b
c
R2
Range
a
b
c
R2
Range
4.031
3.574
2.410
0.898
R = 1.19 – 2.31
1.580
1.846
0.979
0.9626
R = 1.19 – 2.31
Fig. 18. Verification of predictions for the correlation developed for U (i.e. f1). Fig. 15. Variation of liquid fractions for an exemplary case (i.e. Case 17).
Fig. 19. Verification of predictions for the correlation developed for modified L (i.e. f2).
Fig. 16. Liquid fraction variation of U versus the PC model (modified).
797
Solar Energy 188 (2019) 788–798
M.M. Joybari, et al.
correlations (i.e. PC model liquid fraction and storage radius ratio).
energy storage system. Renew. Energy. Katsman, L., Dubovsky, V., Ziskind, G., Letan, R., 2007. Experimental Investigation of Solid-Liquid Phase Change in Cylindrical Geometry. In: ASME/JSME 2007 Thermal Engineering Heat Transfer Summer Conference collocated with the ASME 2007 InterPACK Conference, American Society of Mechanical Engineers, pp. 239–244. Klimeš, L., Mauder, T., Charvát, P., Štětina, J., 2018. Front tracking in modelling of latent heat thermal energy storage: assessment of accuracy and efficiency, benchmarking and GPU-based acceleration. Energy 155, 297–311. Li, C.-Y., Garimella, S.V., Simpson, J.E., 2003. Fixed-grid front-tracking algorithm for solidification problems, part I: method and validation. Numer. Heat Transf., Part B: Fundam. 43 (2), 117–141. Liu, C., Groulx, D., 2014. Experimental study of the phase change heat transfer inside a horizontal cylindrical latent heat energy storage system. Int. J. Therm. Sci. 82, 100–110. Liu, Z., Wang, Z., Ma, C., 2006. An experimental study on the heat transfer characteristics of a heat pipe heat exchanger with latent heat storage Part II: simultaneous charging/ discharging modes. Energy Convers. Manage. 47 (7), 967–991. Liu, Z., Yu, Z.J., Yang, T., Qin, D., Li, S., Zhang, G., Haghighat, F., Joybari, M.M., 2018. A review on macro-encapsulated phase change material for building envelope applications. Build. Environ. Long, J.-Y., Zhu, D.-S., 2008. Numerical and experimental study on heat pump water heater with PCM for thermal storage. Energy Build. 40 (4), 666–672. Mahdi, J.M., Nsofor, E.C., 2018. Solidification enhancement of PCM in a triplex-tube thermal energy storage system with nanoparticles and fins. Appl. Energy 211, 975–986. Mahdi, J.M., Lohrasbi, S., Ganji, D.D., Nsofor, E.C., 2019. Simultaneous energy storage and recovery in the triplex-tube heat exchanger with PCM, copper fins and Al2O3 nanoparticles. Energy Convers. Manage. 180, 949–961. Mat, S., Al-Abidi, A.A., Sopian, K., Sulaiman, M., Mohammad, A.T., 2013. Enhance heat transfer for PCM melting in triplex tube with internal–external fins. Energy Convers. Manage. 74, 223–236. Mirzaei, P.A., Haghighat, F., 2012. Modeling of phase change materials for applications in whole building simulation. Renew. Sustain. Energy Rev. 16 (7), 5355–5362. Murray, R.E., Groulx, D., 2014. Experimental study of the phase change and energy characteristics inside a cylindrical latent heat energy storage system: Part 2 simultaneous charging and discharging. Renew. Energy 63, 724–734. Najafian, A., Haghighat, F., Moreau, A., 2015. Integration of PCM in domestic hot water tanks: optimization for shifting peak demand. Energy Build. 106, 59–64. Omojaro, P., Breitkopf, C., 2014. Investigating and modeling of simultaneous charging and discharging of a PCM heat exchanger. Energy Procedia 48, 413–422. Omojaro, A.P., Breitkopf, C., 2017. Study on solid liquid interface heat transfer of PCM under simultaneous charging and discharging (SCD) in horizontal cylinder annulus. Heat Mass Transf. 1–18. Rathod, M.K., Jyotirmay, B., 2015. Development of correlation for melting time of phase change material in latent heat storage unit. Energy Procedia 75, 2125–2130. Riley, D., Smith, F., Poots, G., 1974. The inward solidification of spheres and circular cylinders. Int. J. Heat Mass Transf. 17 (12), 1507–1516. Seddegh, S., Wang, X., Joybari, M.M., Haghighat, F., 2017a. Investigation of the effect of geometric and operating parameters on thermal behavior of vertical shell-and-tube latent heat energy storage systems. Energy 137, 69–82. Seddegh, S., Wang, X., Henderson, A.D., 2015. Numerical investigation of heat transfer mechanism in a vertical shell and tube latent heat energy storage system. Appl. Therm. Eng. 87, 698–706. Seddegh, S., Wang, X., Henderson, A.D., 2016. A comparative study of thermal behaviour of a horizontal and vertical shell-and-tube energy storage using phase change materials. Appl. Therm. Eng. 93, 348–358. Seddegh, S., Joybari, M.M., Wang, X., Haghighat, F., 2017b. Experimental and numerical characterization of natural convection in a vertical shell-and-tube latent thermal energy storage system. Sustain. Cities Soc. 35, 13–24. Sharifi, N., Faghri, A., Bergman, T.L., Andraka, C.E., 2015. Simulation of heat pipe-assisted latent heat thermal energy storage with simultaneous charging and discharging. Int. J. Heat Mass Transf. 80, 170–179. Solomon, A., 1979. Melt time and heat flux for a simple PCM body. Sol. Energy 22 (3), 251–257. Sparrow, E., Broadbent, J., 1982. Inward melting in a vertical tube which allows free expansion of the phase-change medium. ASME J. Heat Transf. 104 (2), 309–315. Stritih, U., Tyagi, V.V., Stropnik, R., Paksoy, H., Haghighat, F., Mastani Joybari, M., 2018. Integration of passive PCM technologies for net-zero energy buildings. Sustain. Cities Soc. 41, 286–295. Tay, N., Bruno, F., Belusko, M., 2012c. Experimental validation of a CFD and an ε-NTU model for a large tube-in-tank PCM system. Int. J. Heat Mass Transf. 55 (21), 5931–5940. Tay, N., Belusko, M., Bruno, F., 2012a. An effectiveness-NTU technique for characterising tube-in-tank phase change thermal energy storage systems. Appl. Energy 91 (1), 309–319. Tay, N., Belusko, M., Bruno, F., 2012b. Experimental investigation of tubes in a phase change thermal energy storage system. Appl. Energy 90 (1), 288–297. Voller, V.R., Prakash, C., 1987. A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems. Int. J. Heat Mass Transf. 30 (8), 1709–1719. Yang, J., Yang, L., Xu, C., Du, X., 2016. Experimental study on enhancement of thermal energy storage with phase-change material. Appl. Energy 169, 164–176.
• The method can predict the location of the front under SCD outperforming the effective thermal conductivity method. Therefore, the developed method can replace the effective thermal conductivity method especially due to its simplicity. Overall, the users can obtain results within a discrepancy of 15%.
Acknowledgement The authors would like to express their gratitude to Concordia University for the support through Concordia Research Chair – Energy and Environment. References Abdulateef, A.M., Abdulateef, J., Mat, S., Sopian, K., Elhub, B., Mussa, M.A., 2018. Experimental and numerical study of solidifying phase-change material in a triplextube heat exchanger with longitudinal/triangular fins. Int. Commun. Heat Mass Transf. 90, 73–84. Al-abidi, A.A., Mat, S.B., Sopian, K., Sulaiman, M., Mohammed, A.T., 2013a. CFD applications for latent heat thermal energy storage: a review. Renew. Sustain. Energy Rev. 20, 353–363. Al-Abidi, A.A., Mat, S., Sopian, K., Sulaiman, M., Mohammad, A.T., 2013b. Experimental Investigation of Melting in Triplex Tube Thermal Energy Storage. Latest Trends in Renewable Energy and Environmental Informatics. Al-Abidi, A.A., Mat, S., Sopian, K., Sulaiman, M., Mohammad, A.T., 2013c. Experimental study of PCM melting in triplex tube thermal energy storage for liquid desiccant air conditioning system. Energy Build. 60, 270–279. Al-Abidi, A.A., Mat, S., Sopian, K., Sulaiman, M., Mohammad, A.T., 2013d. Internal and external fin heat transfer enhancement technique for latent heat thermal energy storage in triplex tube heat exchangers. Appl. Therm. Eng. 53 (1), 147–156. Al-Abidi, A.A., Mat, S., Sopian, K., Sulaiman, M., Mohammad, A.T., 2013e. Numerical study of PCM solidification in a triplex tube heat exchanger with internal and external fins. Int. J. Heat Mass Transf. 61, 684–695. Al-Abidi, A.A., Mat, S., Sopian, K., Sulaiman, M., Mohammad, A.T., 2014. Experimental study of melting and solidification of PCM in a triplex tube heat exchanger with fins. Energy Build. 68, 33–41. ANSYS Fluent Theory Guide, ANSYS, Inc., Canonsburg, PA, USA, 2013. Bruno, M.B.F., 2008. Design methodology of PCM thermal storage systems with parallel plates. EUROSUN 2008. El-Sawi, A., Haghighat, F., Akbari, H., 2014. Assessing long-term performance of centralized thermal energy storage system. Appl. Therm. Eng. 62 (2), 313–321. Fan, L.-W., Zhu, Z.-Q., Zeng, Y., Lu, Q., Yu, Z.-T., 2014. Heat transfer during melting of graphene-based composite phase change materials heated from below. Int. J. Heat Mass Transf. 79, 94–104. Rubitherm Technologies GmbH, in: Organic PCM, http://www.rubitherm.eu/en/index. php/productcategory/organische-pcm-rt/, Berlin, Germany, 2017. The Copper Tube Handbook, Copper Development Association Inc., New York, NY, USA, 2011. Hasan, A., 1994a. Thermal energy storage system with stearic acid as phase change material. Energy Convers. Manage. 35 (10), 843–856. Hasan, A., 1994b. Phase change material energy storage system employing palmitic acid. Sol. Energy 52 (2), 143–154. Ismail, K.A., Maria das Graças, E., 2003. Melting of PCM around a horizontal cylinder with constant surface temperature. Int. J. Therm. Sci. 42 (12), 1145–1152. Jian-you, L., 2008. Numerical and experimental investigation for heat transfer in triplex concentric tube with phase change material for thermal energy storage. Sol. Energy 82 (11), 977–985. Jones, B.J., Sun, D., Krishnan, S., Garimella, S.V., 2006. Experimental and numerical study of melting in a cylinder. Int. J. Heat Mass Transf. 49 (15), 2724–2738. Joybari, M.M., Haghighat, F., Moffat, J., Sra, P., 2015. Heat and cold storage using phase change materials in domestic refrigeration systems: the state-of-the-art review. Energy Build. 106, 111–124. Joybari, M.M., Haghighat, F., 2016. Natural convection modelling of melting in phase change materials: Development of a novel zonal effectiveness-NTU. 9th International Conference on Indoor Air Quality Ventilation and Energy Conservation In Buildings (IAQVEC 2016), Incheon Songdo, South Korea. Joybari, M.M., Haghighat, F., Seddegh, S., Al-Abidi, A.A., 2017a. Heat transfer enhancement of phase change materials by fins under simultaneous charging and discharging. Energy Convers. Manage. 152, 136–156. Joybari, M.M., Haghighat, F., Seddegh, S., 2017b. Natural convection characterization during melting of phase change materials: development of a simplified front tracking method. Sol. Energy 158, 711–720. Joybari, M.M., Haghighat, F., Seddegh, S., 2017c. Numerical investigation of a triplex tube heat exchanger with phase change material: simultaneous charging and discharging. Energy Build. 139, 426–438. Joybari, M.M., Seddegh, S., Wang, X., Haghighat, F., 2019. Experimental investigation of multiple tube heat transfer enhancement in a vertical cylindrical latent heat thermal
798