Next Article in Journal
Influence of the Size and Shape of Magnetic Core on Thermal Parameters of the Inductor
Next Article in Special Issue
Thermophysical Properties Characterization of Sulphoaluminate Cement Mortars Incorporating Phase Change Material for Thermal Energy Storage
Previous Article in Journal
Pathways for Low-Carbon Transition of the Steel Industry—A Swedish Case Study
Previous Article in Special Issue
Energy Savings in an Office Building with High WWR Using Glazing Systems Combining Thermochromic and Electrochromic Layers
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Thermal Analysis of Pure and Nanoparticle-Enhanced PCM—Application in Concentric Tube Heat Exchanger

School of Mechanical Engineering, National Technical University of Athens, Zografou Campus, 15780 Athens, Greece
*
Author to whom correspondence should be addressed.
Energies 2020, 13(15), 3841; https://doi.org/10.3390/en13153841
Submission received: 8 June 2020 / Revised: 20 July 2020 / Accepted: 24 July 2020 / Published: 27 July 2020

Abstract

:
In this paper, organic phase change materials are modeled and studied both numerically and computationally. The results are compared with those of other research available in the international literature. Regardless of their heat storing capacity, the low heat conductivity of PCMs constitutes a significant obstacle to the further advancement of technologies pertaining to latent heat energy storage, as it hampers the immediate response of systems both at heat storage and at heat recovery. In this work, two types of nanoparticles are used as enhancing media, and the calculations are repeated in terms of melting front and stored energy enhancement. The results show that a 6.5% value for the melt percentage enhancement and a 5.5% value for the heat stored enhancement are exhibited when copper nanoparticles are utilized.

1. Introduction

The constantly increasing emissions of greenhouse gases as well as numerous other gases toxic for human organisms necessitate urgent measures in the direction of energy saving and reducing fossil fuel consumption. Significant contributions to achieving both objectives can be expected by improving the energy performance of buildings to secure lower energy consumption for heating and cooling, and by efficiently utilizing the power generated by means of the thermal storage of untapped heat. In addition, the more efficient utilization of renewable energy sources (RES) shall result in an increase in the share of RES in gross final energy consumption and, consequently, in reduced CO2 emissions.
Among diverse technologies, phase change materials (PCMs) offer considerable advantages, rendering them ideal for thermal energy storage as well as for preserving constant temperatures for longer time spans. PCMs capitalize on the latent heat of fusion and, therefore, store more thermal energy than sensible heat materials. Moreover, during the phase change, PCMs maintain an almost stable temperature, thus contributing to the reduction of the required heat and cooling loads of the buildings [1]. Despite the aforementioned advantages, PCMs require longer time spans to change from the solid to the liquid state and vice versa, due to their low thermal conductivity coefficient. The complexity of the equations describing the phenomenon of the heat flow through the PCM (as a result of the phase change process and the random movement of the fusion-solidification front) does not lend itself to the formation of analytic relations for the calculation of the temperature distribution in the PCM and the heat flow. Consequently, for the purposes of the mathematical analysis and the description of the phenomenon, both the numerical solution and mathematical simulation of the equations are required as well as the experimental calculation of the relevant quantities. Several computational methods have been developed for the accurate simulation of the heat flow through the PCM, describing, with adequate accuracy, the zone where both the liquid and solid phases co-exist (mushy zone). It is typically considered that the material in the mushy zone is solid, partly porous and transforming as the melted material increases. One method of simulation for PCM behavior utilized by several researchers [2,3,4,5] is the enthalpy-porosity method, according to which the enthalpy of the PCM also depends on the ratio of the melted PCM. Another methodology used to solve the heat flow through PCMs by means of finite elements is the lattice Boltzmann method. Both the stability and the efficiency of the method in solving the particular heat transfer problem have been studied by several researchers [6,7,8,9,10]. In addition to the numerical methods used for the solution of the differential equations describing the phenomenon of the heat transfer to and from the PCM, ready-to-use software has also been developed for the simulation of each PCM application. EnergyPlus [11] and TRNSYS [12,13,14] are two commercial software packages widely used for modeling the phase change of the PCM. Further to the numerical solution of the problem arising from each application of PCMs, the calculation of the thermal properties of PCMs (the latent heat of fusion-solidification, specific heat capacity, temperature range of the phase change, and thermal conductivity coefficient) is also imperative as they change during the cooling and heating processes in the material. By utilizing Differential Scanning Calorimetry (DSC), researchers Li et al. [15] calculated the thermal properties of several commercially available PCMs for different values of temperature variation in the chamber. Based on the measurements, the conclusion was that both the temperature at which phase change takes place and the enthalpy variation of the material depend on the rate at which the temperature change occurs. Therefore, further research is required on whether the values of the PCM’s thermal properties match the respective values in each application. Similar were the conclusions reached by researchers M. Iten et al. [16], who studied the effect of specific heat capacity on the phase change process. The results reveal the dependence of the specific heat capacity of the PCM on the rate of temperature change, which may lead to incorrect calculations and thus conclusions if not considered during simulation. Due to the difficulties in the experimental calculation of the thermal properties of PCMs, semi-empirical correlations have been developed to associate the PCM thermal properties with the temperature of the material. G. Ferrer et al. [17] developed empirical correlations to calculate the viscosity and specific heat capacity of fatty acids as well as the viscosity of paraffins [18].
The present paper includes the two-dimensional study of organic PCMs such as paraffin by building a model in Matlab (R2016b, MathWorks) code and by utilizing commercial software, such as Comsol and Fluent (Multiphysics 5.2). The analysis concerns the melting speed and front as well as the heat storage. Moreover, the analysis is repeated with nanoparticles being used as thermal conductivity enhancers so as to accelerate the phase change process.

2. Materials and Methods

2.1. The Examined Domain

The examined domain is depicted in Figure 1. This work regards the analysis of two different geometries. Therefore, Figure 1a shows a rectangular geometry of length 50 mm, which consists of three adiabatic sides and a fourth one subject to a heat source. Figure 1b illustrates the domain used in the concentric tube heat exchanger application in which, apart from the pure PCM analysis, nanoparticles are added in order to increase the low conductivity of the paraffins. The dimensions of the inner and outer tubes are 22 mm and 125 mm, respectively.

2.2. Mathematical Modeling

2.2.1. Pure PCM Mathematical Modeling

The problem of the two-dimensional study, as opposed to the one-dimensional one, requires the solution of all Navier—Stokes equations due to the movement of the liquid PCM caused by the temperature variation in the material. The movement of the PCM facilitates heat flow by natural convection, resulting in the faster melting and distortion of the phase change front. Presented below is the system of equations solving the phase change phenomenon in the material in Cartesian coordinates (Equations (1)–(4)) and cylindrical coordinates (Equations (5)–(8)) [19,20].
ρ t + ( ρ u ) x + ( ρ v ) y = 0
( ρ u ) t + u ( ρ u ) x + v ( ρ u ) y = P x + x ( μ u x ) + y ( μ u y ) S u      
( ρ v ) t + u ( ρ v ) x + v ( ρ v ) y = P y + x ( μ v x ) + y ( μ v y ) S v + ρ g β ( T T m )
( ρ c p T ) t + u ( ρ c p T ) x + v ( ρ c p T ) y = x ( k T x ) + y ( k T y )
ρ t + ( ρ u ) r + ρ u r + ( ρ v ) θ = 0
( ρ u ) t + u ( ρ u ) r + v r ( ρ u ) y v 2 r = P r + 1 r r ( μ r u r ) + 1 r 2 θ ( μ u θ ) μ 2 r 2 v θ μ u r 2 + ρ g β ( T T m ) cos θ S u
( ρ v ) t + u ( ρ v ) r + v r ( ρ v ) y u v r = 1 r P θ + 1 r r ( μ r v r ) + 1 r 2 θ ( μ v θ ) + μ 2 r 2 u θ μ v r 2 + ρ g β ( T T m ) sin θ S v
( ρ c p T ) t + u ( ρ c p T ) r + v r ( ρ c p T ) θ = 1 r r ( r k T r ) + 1 r 2 θ ( k T θ )
In order to zeroize the velocity at the points where the PCM is entirely solid [21], the term S is inserted in the equations for the conservation of momentum, which is proportional to the melted material ratio and given by the following expression:
S = 10 4 · ( 1 e ) 2 0.001 + e 3
Additionally, for the calculation of the pressure field within the liquid PCM, another equation is used [21] to correct the initial hypothetical pressure at any point in time in Cartesian and cylindrical coordinates.
x ( P c o r x S + 1 d t ) + y ( P c o r y S + 1 d t ) = ρ t + ( ρ u ) x + ( ρ v ) y
r ( P c o r r S + 1 d t ) + 1 r ( P c o r r S + 1 d t ) + 1 r θ ( 1 r P c o r θ S + 1 d t ) = ρ t + 1 r ( ρ r u ) r + 1 r ( ρ v ) θ
As illustrated by the equations above, density is calculated using the Boussinesq approximation. This method offers satisfactory results for incompressible flows whose other quantities are stable or temperature-dependent variables. In this particular method, the material density is considered to be stable and a thermal expansion constant is used, β (1/Κ) = 0.001. This method is used to model natural convection provided that the temperature difference is adequately low, i.e.,:
β ( T T ο ) 1
The approximation formula is:
( ρ ρ ο ) g ρ ο   β ( T T ο ) g
The operating density (ρο) is separate for each material, and the operating temperature (Το) is set at 50 °C. In cases when there is great temperature difference or the fluid is compressible, a different correlation between the quantities is required, such as a function correlating temperature with density as a result of experimental data [22].

2.2.2. Effect of Nanoparticle Addition to Phase Change Material

The study of enhanced phase change materials requires the appropriate modification of thermal properties depending on the PCM content in nanoparticles. The density, thermal conductivity, viscosity and heat capacity of the nano-enhanced organic PCM studied above are calculated based on the following expressions [23,24]. The nanoparticles used for analysis are Cu and Al2O3 (Table 1), and the organic PCM is RT50 (Table 2). It should be stressed that for the purposes of the present study, constant nanoparticle properties have been considered while the effect of their size is negligible. Therefore, in the case of nanoparticle-enhanced PCM, the governing equations are the ones described above as well as the following.
k t o t = k p c m k p + 2 k p c m + 2 φ ( k p k p c m ) k p + 2 k p c m φ ( k p k p c m )
ρ t o t = φ ρ p + ( 1 φ ) ρ p c m
μ t o t = 0.983 μ p c m exp ( 12.959 φ p )
c p , t o t = φ ρ p c p , p + ( 1 φ ) ρ p c m c p , p c m ρ t o t
where index p refers to nanoparticles and φ is the volume concentration of the enhanced PCM.

2.3. Solution Procedure of the Phase Change Problem

2.3.1. Numerical Solution Scheme

The method of coupling between the momentum and continuity is of vital importance for both the solution process stability and the total computational time required. As to the problem in question, there was difficulty observed in reaching the continuity equation convergence required, which resulted in either multiple iterations per time step or even divergence among iterations. In general, numerical methods are divided into segregated and coupled methods. In the first type, equations are solved successively, while in the second type, equations are solved simultaneously. The successive solution of equations utilizes some default values or the solution products of the previous cycle and solves each equation in sequence, feeding the results of the solution to the next equation. At the end, each equation is controlled for convergence. In the simultaneous solution of equations, all equations are solved in parallel. The first method has lower memory requirements, since one equation is stored at a time, while for the second method, all equations of the problem need to be registered. However, simultaneous solution tends to be faster. For the solution of the equation system above, a finite-difference time-domain scheme was used, which was solved by means of code developed in Matlab. One basic assumption for the solution of clipped equations is a stable density value for all terms with the exception of the force term in the momentum equations. The solution procedure is illustrated in the following logical diagram (Figure 2).
For the purposes of modeling the phase change phenomenon, a modification of the thermophysical properties of the PCM is required [25,26]. The equations calculating the heat capacity, the density and the thermal conductivity coefficient of the material are given by the equations presented below:
e ( T ) = { 0 , T < T s T T s Δ T , T s < T < T l 1 , T > T l
ρ = ρ s + ( ρ l ρ s ) · e ( T )
k = k s + ( k l k s ) · e ( T )
C p = C p , s + ( c p , l c p , s ) · e ( T ) + L Δ T   · e ( T )
where Δ T is the temperature range of the phase transition and T s ,   T l are the extreme temperatures of the phase transition range.
Finally, Table 3 cites the boundary conditions used in the solution in Cartesian and cylindrical coordinates.

2.3.2. Computational Solution Scheme

Further to the numerical solution of the Navier—Stokes equations, which describes the phase change phenomenon, simulations of the phase change phenomenon were also performed by means of the COMSOL and Fluent software. For the purposes of the solution for the phenomenon, appropriate modifications to the thermal properties were made to separate the liquid from the solid phases [26,27].
μ = μ l ( 1 + A ( T ) )
A ( T ) = 10 5 ( 1 e ( T ) ) 2 ( e ( T ) 3 + 10 3 )
L · D ( T ) = L · e ( T T p c m ) 2 Δ T 2 π · Δ T 2       , T s < T < T l
As is generally required in computational methods, within the surface or domain to be studied, computational nodes are specified by means of various methods, thus determining a finite number of computational cells, for each of which a series of calculations is performed. The formation of the appropriate mesh is particularly significant for the proper simulation of the behavior of the various models and phenomena to be studied in each case. In general, the denser the mesh is, the better the approximation but with higher demands on computational resources (CPU time and memory). Therefore, standard practice provides for the determination of meshes that are denser near the critical regions (complex geometries, orifices, boundary layers and plumes) and thinner in areas expected to be of less interest. In the particular case studied, the natural convection occurring due to the temperature difference in the paraffin melt highlights the importance of the mesh for the results. The mesh was defined in terms of triangle cells and was structured using the patch-dependent method. According to this method, the mesh is determined by considering the number of nodes defined at the edge-curves and the maximum and minimum distances defined between the nodes [28].

3. Results and Discussion

The PCM selected for the two-dimensional analysis of the phase change was paraffin RT50, which was also used for the one-dimensional analysis. A constant time step at 0.1 s was defined for the solution of the problem. For the completion of each time step and the attainment of convergence, constantly increasing iterations are required. This is due to the fact that as PCM melting proceeds, more and more computational nodes or cells (in the case of COMSOL) are affected by natural convection and acquire the speed and properties of fluid, thus leading to an increase in residuals and, consequently, in the number of iterations required. Therefore, the computational time corresponding to each real second increases and the model solution decelerates. It is worth mentioning characteristically that while, initially, the iterations required are 20–30 per time step, after 900 s and 9000 steps, the iterations required are over 100 per iteration. The total time required for each of the cases so as to reach 30 min in the evolution of the phenomenon is over 5 days. Moreover, in the case of the pure PCM, in which the kinematic connectivity is lower and higher speeds are developed, the convergence of the continuity equation is even harder to attain and the total iterations per time step further increase.

3.1. Grid Independence

A solution by a numerical or computational method is controlled in terms of validity by the grid independence in the computational domain. In the case of the solution code in Matlab, the temperature variation at one particular position within the material is used for selecting the density of the grid. Thus, in the case of the numerical solution, an 81 × 81 grid is selected as shown in Figure 3.
As far as the computational solution by means of software is concerned, the effective simulation of models in CFD software (Multiphysics 5.2) is largely dependent on the mesh density and the time step. Further to setting these parameters with a view to calculation convergence, it is also imperative that the grid density and the size of the time step do not affect, within certain boundaries, the results of the model solution. For this particular reason, prior to the final solution of the model, through trials, with increasingly denser grids, the number of nodes is defined, a further increase in which does not significantly affect the results [28]. Figure 4 presents the melt percentage as function of time for the different grids as well as the evolution over time of the percentage variation from the melt percentage with the highest density. The grids studied exhibit densities of 70,000 and 100,000 computational nodes, respectively. It can be observed that out of the two grids, the one with the fewest nodes exhibits a significant difference in relation to the dense grid, which actually increases over time. On the other hand, the 100,000-node grid exhibits less difference in relation to the dense grid, which actually remains relatively stable over time. Based on the above, the 100,000-node grid was selected as being both accurate and computationally cheaper.
Figure 5 presents the evolution of the melt percentage over time for different time steps and their difference from the smallest time step controlled, which was 0.01 s. The time step of 0.01 s is the least temporal step that can be controlled in such problems; thus, it is used as a reference in our approach. Through this step, Figure 5 was created, as it depicts the deviation of the solution from the reference one for other time steps. The discrepancy is minimized for time step 0.1 s, which was finally selected. As for the selection of the grid and the time step, the cost in terms of computational resources is also considered.

3.2. Rectangular Cross Section Analysis

The geometry selected for the accuracy of the computational model is rectangular, which is the most common geometry for PCM study. The dimensions of the cross section are 5 cm × 5 cm, while the field depth is not defined. The validation of the COMSOL model is accomplished by means of comparison with the available literature. The simulation parameters are adopted by the respective work in Souyafane et al. [29], and the COMSOL model in which the modified properties are inserted runs until 5000 s. The complexity of the problem is illustrated by the fact that in order for the model to store as a result only the melt front up to the time of 5000 s, 98 h of continuous simulation was required. Figure 6 represents the melting front as isothermal curves and includes the comparison with the literature results. As observed, there is agreement between the two methods, which renders the modified COMSOL model reliable.
Following the validation of the COMSOL model, the simulation parameters are inserted in the Matlab mathematical model. The calculations cover the first 250 s and have a total duration of 16.8 h, which is indicative of the problem complexity. Figure 7 presents the melt front for material RT50 of the geometry selected. The buoyancy effect is obvious and deforms the melt front, creating an area in the northwestern corner, which moves at higher velocity than the rest of the front. The melting front is in complete agreement with the available literature, which renders the solution of the Matlab developed model acceptable. In the case in which the effect of gravity through the buoyancy flow is not taken into consideration, the front moves in uniform flow and the one-dimensional model would suffice for the analysis of the phase change materials. Buoyancy creates a velocity field, the maximum value of which is 10−3 m/s. In order to illustrate the evolution of the phenomenon, the mathematical model is again used up to the time of 7200 s with a thinner grid so as to slightly reduce the time required to reach the solution. Figure 7b represents the melt front two hours after the onset of the process. Due to natural convection, the melt PCM flows at the upper side of the cross section, thus accelerating non-uniform melting. The total duration of the simulation for the material studied up to the specific time is 160.6 h.
The validity of the developed numerical model is confirmed through comparison with the respective results calculated by the COMSOL model in Figure 8. The deviation observed in the first seconds and for the low values of the melt percentage is due to the density of the grid selected for the numerical scheme in Matlab. The rapid front expansion for the areas near the heat source can be easily observed. Then, the effect of the heat source diminishes and the front evolves at a low rate. For the same time span, the heat stored in the PCM is also represented. The heat is determined by the PCM surface integral, and the results are represented by field depth. The relatively low values of heat stored as shown in Figure 9 are justified by the fact that they correspond to the heat charging of the material with sensible heat. The discrepancy between the numerical solution and the COMSOL solution is due to the reasons already stated for the case of the melt percentage.

3.3. Concentric Tube Heat Exchange Analysis

3.3.1. Pure PCM Results

The geometry studied consists of a concentric tube heat exchanger with tube diameters of Φ125 and Φ22, respectively. The phase change material RT50 is located in between the concentric tubes while the high temperature (80 °C) is applied to the internal tube. The analysis was conducted in the Fluent environment with the modification of properties similar to those effectuated in COMSOL. The following figures represent the temporal development of the phenomenon and, in particular, the evolution of the melt over the thermal field and the velocity fields for different points in time. Figure 10 presents the variation in the paraffin temperature up to 800 s, while Figure 11 presents the isothermal curves for the whole domain at 900 s. Up to the first 100 s, conduction is the driving force; therefore, the isothermal curves are uniform. After that, gravitational forces and buoyance are starting to affect the phase change and heat is transferred from the source to the material, not only through conduction but via convection as well. Small vortices are created, which enhance the natural convection phenomenon. The shape of the isothermal curve is also typical and refers to the shape of the melt as a result of the buoyancy and the vortices produced by the velocity field.
Farhadi et al. [30] studied the melting of organic PCM with a two-dimensional model of the same geometry as the one analyzed here. The melting formulation of this work is in agreement with the results of Farhadi.
The following figures represent the temporal evolution of the velocity field as well as magnified areas from the field at the point 900 s. Figure 12 illustrates the mode in which velocities are developed as a result of natural convection in the melt as well as the significant differentiation of the measured velocity depending on the area and given the specific scale. In every case, the measured velocity acquires low values—typical of the buoyancy flow—which vary from 10−5 m/s in the first 100 s up to 10−3 m/s at 900 s. It is also important to stress the observed symmetry as to the vertical axis y in Figure 13.
Within the areas at which melt circulation occurs, there are areas with near zero velocity at specific points. This can be observed at Figure 14, which represents a magnified area of the vector velocity field at point t = 900 s; it illustrates the vortex formed due to the natural circulation of the melt, which justifies the zero velocity at the center. As a result of this circulation, in which the hot and light part of the melt rises until it reaches the solid edge or loses its heat load and cools down, several similar vortices and equal zero velocity areas are developed. Finally, the highest velocities per meter are observed at the confluence of individual vortices.
Finally, Figure 15 represents the temporal variation in the energy stored in the mass of the phase change material. This heat results from the completion of the enthalpy, the melt percentage fluctuations at the incremental surface of the material studied, and the ensuing conversion to the total surface of the PCM.

3.3.2. Nanoparticle-Enhanced PCM

The following figures represent the effects of adding nanoparticles on the melt percentage (Figure 16) as well as the heat stored in the enhanced PCM (Figure 17). The percentage increase in the melt appears to be unstable at the beginning of the process because the melt quantity in the first seconds is infinitely small, which causes abrupt fluctuation. With the increase in the melt quantity comes the stabilization of the percentage or at least a smoother variation. When melting begins and heat transfer by conduction prevails, the enhanced materials with the highest thermal conductivity, i.e., PCM-Cu (3%) and PCM-Al2O3 (3%), have better performance in terms of both heat storage and material melting. After 300 s of operation and the production of an efficient melt quantity, natural circulation intensifies and natural conduction becomes the primary heat transfer mechanism. The outcome is that copper nanoparticles are mostly advantageous in enhancing the thermal behavior of the PCM, whereas PCM-Al2O3 (3%), from the establishment of natural convection and onwards, falls behind in terms of melt growth even when compared with PCM-Al2O3 (1%). The fact that PCM-Al2O3 is less effective than PCMs-Cu may be attributed mainly to its lower thermal conductivity. One major parameter influencing the speed of melting and heat storage is the kinematic viscosity of the enhanced material, which increases with the addition of nanoparticles. The increase in kinematic viscosity inhibits the growth of melt circulation, thus suspending natural convection and, consequently, heat transfer from the hot wall to the enhanced PCM. Therefore, heat is transferred faster in the lower heat conductivity PCM-Cu (1%) than in PCM- Al2O3 (3%), which is attributed to its low kinematic viscosity, especially after 300 s, when the circulation of the superheated melt has reached satisfactory levels. The importance of kinematic viscosity has even been stressed in the literature but is also prominent in the slight prevalence of PCM-Al2O3 (1%) over PCM-Al2O3 (3%).

4. Conclusions

Considering the two-dimensional analysis accomplished by means of developed code and available software, and which focuses on both pure organic PCM and the enhanced material with metal and metal oxide nanoparticles, the main conclusions are summarized as follows:
  • The complexity of the problem, especially in the “mushy zone” area, necessitates particularly long simulation periods and increased computing power for the solution. For example, the evolution of the phenomenon up to 7200 s requires approximately 160.6 h of simulation.
  • Buoyancy induces velocity fields with measured velocity of about 10−3 m/s. This natural circulation of the melt intensifies the heat transfer and the phase change as it transfers heat from the melt to the solid by a combination of natural convection and conduction.
  • The velocity field takes the form of vortex with areas of zero velocity (vortex center) and areas of maximum velocity (confluence with neighboring vortex).
  • The addition of Cu and Al2O3 nanoparticles to the compositions studied enhances both the melt percentage and the heat stored in the material. However, this effect declines over time for the case of the enhanced PCM bearing Al2O3 nanoparticles. On the other hand, for the case of copper nanoparticles, this enhancement exhibits intense and irregular fluctuations for the first seconds of the process and then maintains stable values of about 6.5% for the melt percentage and 5.5% for the heat stored.
  • The quantity of the nanoparticles added cannot be increased infinitely, and the reason for this is the kinematic viscosity of the enhanced material, which increases with the addition of nanoparticles. The increase in kinematic viscosity inhibits the growth of melt circulation, thus suspending natural convection and, consequently, heat transfer from the hot wall to the enhanced PCM.

Author Contributions

Conceptualization, M.T.N. and I.P.K.; Methodology, M.T.N.; Software, M.T.N.; Validation, M.T.N. and I.P.K.; Formal Analysis, M.T.N.; Investigation, M.T.N. and I.P.K. Resources, M.T.N.; Data Curation, M.T.N. and I.P.K.; Writing-Original Draft Preparation, M.T.N. and I.P.K.; Writing-Review & Editing, I.P.K.; Visualization, M.T.N.; Supervision, I.P.K.; Project Administration, I.P.K.; Funding Acquisition, M.T.N. All authors have read and agreed to the published version of the manuscript.

Funding

This research is co-financed by Greece and the European Union (European Social Fund—ESF) through the Operational Programme «Human Resources Development, Education and Lifelong Learning» in the context of the project “Reinforcement of Postdoctoral Researchers—2nd Cycle” (MIS-5033021), implemented by the State Scholarships Foundation (ΙΚΥ).
Energies 13 03841 i001

Acknowledgments

Michael Nitsas would like to thank the State Scholarships Foundation (IKY) for the financial support (MIS-5033021). Both authors acknowledge with appreciation the Laboratory of Applied Thermodynamics, Thermal Engineering Section, School of Mechanical Engineering, National Technical University of Athens, for the support of the work on which this paper is based.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

Latin letters
u, v[m/s]Velocity components
g[m/s2]Gravitational acceleration
D, d[m]Diameter
Τ[K], [οC]Temperature
k[W/mK]Thermal conductivity
Cp[J/kgK]Specific heat
p[Pa]Pressure
A[m2]Surface
q[W/m2]Thermal flux
LmLength
h[W/m2K]Convection coefficient
X, YmCoordinates
m[kg/s]Mass flow rate
e[-]Melting fraction
E[-]Enhancement
r[m]Radius
h[KJ/kg]Specific enthalpy
L[KJ/kg]Latent heat of melting
S[m]Displacement of melting front
Q[J], [W]Heat
Greek
v [m2/s]Kinematic viscosity, v = μ/ρ
βν[1/Κ]Volume expansion coefficient
ρ[kg/m3]Density
φ [-]Volumetric nanoparticle concentration
μ[Pas]Dynamic viscosity
Subscripts and Superscripts
p Nanoparticles
s Solid
l Liquid
tot Total
melt Melting
pcm Phase change
stored Stored
comp Composite
cPCM composite PCM
hybrid Hybrid
Abbreviations
PCM Phase change material

References

  1. Fleischer, A.S. Thermal Energy Storage Using Phase Change Materials; Springer: New York, NY, USA, 2015; ISBN 978-3-319-20922-7. [Google Scholar]
  2. Chakraborty, P.R. Enthalpy porosity model for melting and solidification of pure-substances with large difference in phase specific heats. Int. Commun. Heat Mass Transf. 2017, 81, 183–189. [Google Scholar] [CrossRef]
  3. Voller, V.; Prakash, C. A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems. Int. J. Heat Mass Transf. 1987, 30, 1709–1719. [Google Scholar] [CrossRef]
  4. Niezgoda-Żelasko, B. The Enthalpy-porosity Method Applied to the Modelling of the Ice Slurry Melting Process During Tube Flow. Procedia Eng. 2016, 157, 114–121. [Google Scholar] [CrossRef] [Green Version]
  5. Brent, A.D.; Voller, V.R.; Reid, K.J. Enthalpy-Porosity technique for modeling convection-diffusion phase change: Application to the melting of a pure metal. Numer. Heat Transf. 1988, 13, 297–318. [Google Scholar] [CrossRef]
  6. Taghilou, M.; Talati, F. Numerical investigation on the natural convection effects in the melting process of PCM in a finned container using lattice Boltzmann method. Int. J. Refrig. 2016, 70, 157–170. [Google Scholar] [CrossRef]
  7. Talati, F.; Taghilou, M. Lattice Boltzmann application on the PCM solidification within a rectangular finned container. Appl. Therm. Eng. 2015, 83, 108–120. [Google Scholar] [CrossRef]
  8. Siebert, D.N.; Hegele, L.; Philippi, P. Lattice Boltzmann equation linear stability analysis: Thermal and athermal models. Phys. Rev. E 2008, 77, 026707. [Google Scholar] [CrossRef]
  9. Yuan, P. Thermal Lattice Boltzmann Two Phase Flow Model for Fluid Dynamics. Ph.D. Thesis, University of Pittsburgh, Pittsburgh, PA, USA, 2005. [Google Scholar]
  10. Bawazeer, S. Stability and Accuracy of Lattice Boltzmann Method. Master’s Thesis, University of Calgary, Calgary, AB, Canada, 2013. [Google Scholar]
  11. Tabares-Velasco, P.C.; Christensen, C.; Bianchi, M. Verification and validation of EnergyPlus phase change material model for opaque wall assemblies. Build. Environ. 2012, 54, 186–196. [Google Scholar] [CrossRef]
  12. Kuznik, F.; Virgone, J.; Johannes, K. Development and validation of a new TRNSYS type for the simulation of external building walls containing PCM. Energy Build. 2010, 42, 1004–1009. [Google Scholar] [CrossRef] [Green Version]
  13. Al-Saadi, S.N.; Zhai, Z.Q. A new validated TRNSYS module for simulating latent heat storage walls. Energy Build. 2015, 109, 274–290. [Google Scholar] [CrossRef] [Green Version]
  14. Lu, S.; Liu, S.; Huang, J.; Kong, X. Establishment and experimental verification of PCM room’s TRNSYS heat transfer model based on latent heat utilization ratio. Energy Build. 2014, 84, 287–298. [Google Scholar] [CrossRef]
  15. Li, L.; Yu, H.; Wang, X.; Zheng, S. Thermal analysis of melting and freezing processes of phase change materials (PCMs) based on dynamic DSC test. Energy Build. 2016, 130, 388–396. [Google Scholar] [CrossRef]
  16. Iten, M.; Liu, S.; Shukla, A.; Da Silva, P.D. Investigating the impact of Cp-T values determined by DSC on the PCM-CFD model. Appl. Therm. Eng. 2017, 117, 65–75. [Google Scholar] [CrossRef]
  17. Ferrer, G.; Barreneche, C.; Palacios, A.; Solé, A.; Fernández, A.I.; Cabeza, L.F. Empirical equations for viscosity and specific heat capacity determination of fatty acids. J. Energy Storage 2017, 10, 20–27. [Google Scholar] [CrossRef] [Green Version]
  18. Ferrer, G.; Gschwander, S.; Solé, A.; Barreneche, C.; Fernández, A.I.; Schossig, P.; Cabeza, L.F. Empirical equation to estimate viscosity of paraffin. J. Energy Storage 2017, 11, 154–161. [Google Scholar] [CrossRef] [Green Version]
  19. Lopez, J.; Marques, F.; Shen, J. An Efficient Spectral-Projection Method for the Navier–Stokes Equations in Cylindrical Geometries. J. Comput. Phys. 2002, 176, 384–401. [Google Scholar] [CrossRef] [Green Version]
  20. Ahmad, R.A. Steady-State Numerical Solution of the Navier-Stokes and Energy Equations around a Horizontal Cylinder at Moderate Reynolds Numbers from 100 to 500. Heat Transf. Eng. 1996, 17, 31–81. [Google Scholar] [CrossRef]
  21. Galione, P.A.; Lehmkuhl, O.; Rigola, J.; Oliva, A. Fixed-Grid Modeling of Solid-Liquid Phase Change in Unstructured Meshes Using Explicit Time Schemes. Numer. Heat Transf. Part B Fundam. 2013, 65, 27–52. [Google Scholar] [CrossRef]
  22. Shmueli, H.; Ziskind, G.; Letan, R. Melting in a vertical cylindrical tube: Numerical investigation and comparison with experiments. Int. J. Heat Mass Transf. 2010, 53, 4082–4091. [Google Scholar] [CrossRef]
  23. Yang, L.; Xu, X. A renovated Hamilton–Crosser model for the effective thermal conductivity of CNTs nanofluids. Int. Commun. Heat Mass Transf. 2017, 81, 42–50. [Google Scholar] [CrossRef]
  24. Kavusi, H.; Toghraie, D. A comprehensive study of the performance of a heat pipe by using of various nanofluids. Adv. Powder Technol. 2017, 28, 3074–3084. [Google Scholar] [CrossRef]
  25. Lamberg, P.; Lehtiniemi, R.; Henell, A.-M. Numerical and experimental investigation of melting and freezing processes in phase change material storage. Int. J. Therm. Sci. 2004, 43, 277–287. [Google Scholar] [CrossRef]
  26. Murray, R.; Groulx, D. Modeling Convection during Melting of a Phase Change Material. In Proceedings of the COMSOL Conference, Boston, MA, USA, 13–15 October 2011. [Google Scholar]
  27. Samara, F.; Groulx, D.; Biwole, P. Natural Convection Driven Melting of Phase Change Material: Comparison of Two Methods. In Proceedings of the COMSOL Conference, Boston, MA, USA, 3–5 October 2012. [Google Scholar]
  28. ANSYS. Ansys Fluent 12.0 User’s Guide; ANSYS: Pittsburgh, PA, USA, 2009. [Google Scholar]
  29. Souayfane, F.; Biwole, P.H.; Fardoun, F. Melting of a phase change material in presence of natural convection and radiation: A simplified model. Appl. Therm. Eng. 2018, 130, 660–671. [Google Scholar] [CrossRef]
  30. Darzi, A.R.; Farhadi, M.; Sedighi, K. Numerical study of melting inside concentric and eccentric horizontal annulus. Appl. Math. Model. 2012, 36, 4080–4086. [Google Scholar] [CrossRef]
Figure 1. The analysis domain: (a) Rectangular domain of L = 5 cm with three adiabatic walls; (b) Concentric tube heat exchanger with inner tube diameter d = 22 mm and outer tube diameter D = 125 mm.
Figure 1. The analysis domain: (a) Rectangular domain of L = 5 cm with three adiabatic walls; (b) Concentric tube heat exchanger with inner tube diameter d = 22 mm and outer tube diameter D = 125 mm.
Energies 13 03841 g001
Figure 2. Flow diagram of the numerical scheme.
Figure 2. Flow diagram of the numerical scheme.
Energies 13 03841 g002
Figure 3. Grid density of the domain used in the numerical solution.
Figure 3. Grid density of the domain used in the numerical solution.
Energies 13 03841 g003
Figure 4. Melting fraction variation as a function of time for the examined node densities used in the computational solution.
Figure 4. Melting fraction variation as a function of time for the examined node densities used in the computational solution.
Energies 13 03841 g004
Figure 5. Melting fraction variation for various time steps.
Figure 5. Melting fraction variation for various time steps.
Energies 13 03841 g005
Figure 6. COMSOL results (a) and comparison with available data from literature (b) [29].
Figure 6. COMSOL results (a) and comparison with available data from literature (b) [29].
Energies 13 03841 g006
Figure 7. Melting front of RT50, (a) for the first 250 s and (b) at t = 7200 s.
Figure 7. Melting front of RT50, (a) for the first 250 s and (b) at t = 7200 s.
Energies 13 03841 g007
Figure 8. Melting fraction evolution for RT50.
Figure 8. Melting fraction evolution for RT50.
Energies 13 03841 g008
Figure 9. Stored energy during the interval 0–250 s.
Figure 9. Stored energy during the interval 0–250 s.
Energies 13 03841 g009
Figure 10. Variation of temperature field for 100 s (a), 300 s (b), 500 s (c) and 800 s (d).
Figure 10. Variation of temperature field for 100 s (a), 300 s (b), 500 s (c) and 800 s (d).
Energies 13 03841 g010
Figure 11. Temperature field at 900 s.
Figure 11. Temperature field at 900 s.
Energies 13 03841 g011
Figure 12. Velocity field at 100 s (a), 300 s (b), 500 s (c) and 700 s (d).
Figure 12. Velocity field at 100 s (a), 300 s (b), 500 s (c) and 700 s (d).
Energies 13 03841 g012
Figure 13. Velocity field at t = 900 s.
Figure 13. Velocity field at t = 900 s.
Energies 13 03841 g013
Figure 14. Magnification of Figure 16—formation of molten material vortices.
Figure 14. Magnification of Figure 16—formation of molten material vortices.
Energies 13 03841 g014
Figure 15. Stored energy variation up to t = 900 s.
Figure 15. Stored energy variation up to t = 900 s.
Energies 13 03841 g015
Figure 16. Melting front enhancement with the nanoparticles’ utilization.
Figure 16. Melting front enhancement with the nanoparticles’ utilization.
Energies 13 03841 g016
Figure 17. Thermal energy storage enhancement with the nanoparticles’ utilization.
Figure 17. Thermal energy storage enhancement with the nanoparticles’ utilization.
Energies 13 03841 g017
Table 1. Nanoparticle properties.
Table 1. Nanoparticle properties.
Nanoparticles ρ   ( k g m 3 ) k   ( W m K ) C p ( J k g K )
Al 2 O 3 360036765
Cu 8933400385
Table 2. Simulation parameters.
Table 2. Simulation parameters.
PropertyRT50
Phase change temperature T p c m = 50 °C
Δ T = 45 51 °C
Ts45 °C
Latent heat L = 125   kJ kg
Density (solid state) ρ s = 880   kg m 3
Density (liquid state) ρ l = 760   kg m 3
Specific heat (both phases) c p = 2.0   kJ kgK
Conductivity k = 0.2   W mK
Table 3. Initial and boundary conditions.
Table 3. Initial and boundary conditions.
Cartesian CoordinatesCylindrical Coordinates
T i n = 1 = T s T i n = 1 = T s
T i = N I n = T s T r i = N I n = 0
T i = 1 n = 343   K T i = 1 n = 353   K
T x j = 1 , N J = 0 T r j = 1 , N J = 0
u w a l l = 0 u w a l l = 0

Share and Cite

MDPI and ACS Style

T. Nitsas, M.; P. Koronaki, I. Thermal Analysis of Pure and Nanoparticle-Enhanced PCM—Application in Concentric Tube Heat Exchanger. Energies 2020, 13, 3841. https://doi.org/10.3390/en13153841

AMA Style

T. Nitsas M, P. Koronaki I. Thermal Analysis of Pure and Nanoparticle-Enhanced PCM—Application in Concentric Tube Heat Exchanger. Energies. 2020; 13(15):3841. https://doi.org/10.3390/en13153841

Chicago/Turabian Style

T. Nitsas, M., and I. P. Koronaki. 2020. "Thermal Analysis of Pure and Nanoparticle-Enhanced PCM—Application in Concentric Tube Heat Exchanger" Energies 13, no. 15: 3841. https://doi.org/10.3390/en13153841

APA Style

T. Nitsas, M., & P. Koronaki, I. (2020). Thermal Analysis of Pure and Nanoparticle-Enhanced PCM—Application in Concentric Tube Heat Exchanger. Energies, 13(15), 3841. https://doi.org/10.3390/en13153841

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop