Thermal Energy Storage and Heat Transfer of Nano-Enhanced Phase Change Material (NePCM) in a Shell and Tube Thermal Energy Storage (TES) Unit with a Partial Layer of Eccentric Copper Foam

Thermal energy storage units conventionally have the drawback of slow charging response. Thus, heat transfer enhancement techniques are required to reduce charging time. Using nanoadditives is a promising approach to enhance the heat transfer and energy storage response time of materials that store heat by undergoing a reversible phase change, so-called phase change materials. In the present study, a combination of such materials enhanced with the addition of nanometer-scale graphene oxide particles (called nano-enhanced phase change materials) and a layer of a copper foam is proposed to improve the thermal performance of a shell-and-tube latent heat thermal energy storage (LHTES) unit filled with capric acid. Both graphene oxide and copper nanoparticles were tested as the nanometer-scale additives. A geometrically nonuniform layer of copper foam was placed over the hot tube inside the unit. The metal foam layer can improve heat transfer with an increase of the composite thermal conductivity. However, it suppressed the natural convection flows and could reduce heat transfer in the molten regions. Thus, a metal foam layer with a nonuniform shape can maximize thermal conductivity in conduction-dominant regions and minimize its adverse impacts on natural convection flows. The heat transfer was modeled using partial differential equations for conservations of momentum and heat. The finite element method was used to solve the partial differential equations. A backward differential formula was used to control the accuracy and convergence of the solution automatically. Mesh adaptation was applied to increase the mesh resolution at the interface between phases and improve the quality and stability of the solution. The impact of the eccentricity and porosity of the metal foam layer and the volume fraction of nanoparticles on the energy storage and the thermal performance of the LHTES unit was addressed. The layer of the metal foam notably improves the response time of the LHTES unit, and a 10% eccentricity of the porous layer toward the bottom improved the response time of the LHTES unit by 50%. The presence of nanoadditives could reduce the response time (melting time) of the LHTES unit by 12%, and copper nanoparticles were slightly better than graphene oxide particles in terms of heat transfer enhancement. The design parameters of the eccentricity, porosity, and volume fraction of nanoparticles had minimal impact on the thermal energy storage capacity of the LHTES unit, while their impact on the melting time (response time) was significant. Thus, a combination of the enhancement method could practically reduce the thermal charging time of an LHTES unit without a significant increase in its size.


Introduction
Eighty-one percent of energy consumption was produced from fossil fuel resources in 2017 [1]. Energy lost as waste heat accounts for 12.2% of global energy consumption [2,3]. Renewable energy is under rapid development, so it can replace fossil fuel energy and reduce global warming gas emissions. However, different types of renewable energy are generally intermittent and fluctuant, so it is impossible to use them continuously. Furthermore, there is a discrepancy between the supply and the demand which must be reconciled. One potential solution is to store surplus energy as heat and release it when there is a deficit in power generation. This is known as thermal energy storage (TES). One way to achieve this is to use excess energy to force a material to undergo an endothermic phase change, which can be reversed to release latent heat when needed. This is the concept of latent heat thermal energy storage (LHTES).
One of the advantages of LHTES is its large capacity to store energy, compared with alternative concepts such as sensible heat TES. However, the concept suffers from a low thermal response, which is the main disadvantage in commercialization and large-scale applications of LHTES. The heat transfer efficiency is limited, and the charge/discharge process is lengthened in most materials because of poor thermal conductivity.
To enhance the thermal response of LHTES, various techniques have been proposed, such as the use of a finned unit to increase total surface areas [4,5], direct contact heat transfer [6], and the modification of materials with a porous matrix [7,8]. Several studies have demonstrated that the time required to store and discharge heat is reduced by using foams, such as a copper foam bonded to a material that undergoes the phase change [9], a paraffin/copper foam and paraffin/nickel foam composite [10], and conductive foams and finned pipes [11].
During charging, the primary heat transfer mechanisms are conduction and (natural) convection [12]. Therefore, the performance of any system is dependent on the geometry. For example, in a shell-and-tube heat exchanger, in which a central tube extracts a head from a molten material held in a container ("shell"), more heat is transferred when the axis of both the shell and the tube is horizontal, but the maximum rate of heat transfer is highest when both were vertical [13]. Another study showed that if enough molten material is present in a horizontal shell with a finned copper tube, convection dominates heat transfer [14].
A very important parameter in the design of shell-and-tube TES systems is the eccentricity of the tube [15]. In a system that uses a horizontal shell of a storage material, a reduction in eccentric distance leads to an improvement in the melting rate and an increase in the rate of heat transfer by natural convection [16]. The orientation of the tube can be changed independently of the shell: the orientation of the tube in a horizontal shell system impacts the rate at which heat can be stored (the charging speed) [17].
Using an extended surface or fins can also improve the thermal performance of LHTES systems by providing a higher surface area across which heat may be transferred. Numerous researchers have studied potential configurations of finned heat exchangers such as tree-shaped fins [18], pin fin heat sinks [19,20], and plate-finned heat exchangers [21].
A plate-fin unit with a length of 2 mm for rapid heat storage/release using paraffin was studied numerically. When temperature differences in the system are less than 20 • C, the fins are vital to the energy storage performance [22]. In a fin-and-tube heat exchanger, the use of fins decreases the time for both melting and the solidification of the storage material, independent of the flow regime and pitch of the fins [23]. In a plate-fin-type system, it has been found that the Stefan number, the Rayleigh number, the Nusselt number, and the Fourier number are important factors [24].
Heat sinks with pin fins and materials that undergo phase changes have been utilized widely in the cooling of electronic products. The base temperature and the heat exchanger operating time can be decreased significantly by raising the number, the height, and the thickness of the fins [25,26].
Metal foams could increase the thermal conductivity in LHTES and thereby enhance heat transfer. In one study using expanded graphite and metal foams, an interconnected structure of metal foams enhances heat transfer but reduces convection in the molten material [27]. Increasing the porosity linearly from the bottom of the system to the top can shorten the time required to melt the storage material and improve the heat transfer performance, compared to the case of a constant porosity of the investigated melting in a phase change material (PCM) saturated in an irregular geometry. Mesalhy et al. [28] investigated an irregular geometry filled with a high thermal-conductive porous matrix. The authors showed an increase in the melting rate when the porosity of the matrix decreases, but the convection is hindered. Overall, it is concluded from the literature that the use of a solid matrix with high thermal conductivity and high porosity is the best technique to improve storage response in LHTES.
The net useful energy during charging and discharging periods has been estimated by performing energy and exergy analyses of shell-and tube-type LHTES units for a medium temperature solar thermal power plant (~200 • C) [29]. Similarly, the behavior of a LHTES system in the presence and in the absence of an aluminum foam has been studied. It was shown that the presence of foam accelerates the phase change and the heat transfer of the LHTES system and that the melting time is shown to increase when the porosity of the foam increases [30].
Adding nanoparticles to materials that undergo a phase change in a porous energy storage system decreases the solidification time by up to 23.5%. However, the solidification time decreases by 14% when adding nanoparticles with volume fractions up to 4% [31]. For example, by adding boron nitride nanoparticles to calcium chloride hexahydrate, a thermal conductivity can be 71.9% higher than that with calcium chloride hexahydrat alone. However, the specific heat capacity and the latent heat of fusion decrease 60.9% and 11.1%, respectively, when the nanoparticles are added. Incorporating alumina nanoparticles to a material inside a triplex tube has little effect at the beginning of the solidification process. However, the rate of solidification increases with the volume fraction of nanoparticles [32].
The use of a metal foam, fins, and nanoparticles were each studied in the thermal management of a lithium-ion battery. The foam gave a bigger performance boost than either the fins or the nanoparticles [33]. Separately, a multiple-segment metal foam with an average porosity of 0.95 and an incorporation of 5 vol.% nanoparticles in a tube-andshell LHTES system reduces solidification time by up to 94%, compared to the case of a single pure storage material. The effect is dependent on the number of cascaded foam segments and the number of different materials that are able to undergo a phase change [34]. In a separate study, the simultaneous use of a porous foam and nanoparticles reduces solidification time by 21.4%, compared to the use of a pure material [35].
Metal foams can improve the response time of LHTES units. Using nanoparticles can also improve thermal conductivity and heat transfer. However, a porous matrix limits natural convection flows. A dispersion of nanoparticles raises the dynamic viscosity [36] and, consequently, reduces natural convection flows. Moreover, the addition of metal foams and nanoadditives will increase the masses of LHTES units. In addition, since they do not contribute to latent heat energy storage, the overall heat capacity of the LHTES units can be reduced. The proper design of an LHTES system is a crucial task to benefit from the heat transfer improvement of metal foams and nanoadditives but to avoid overweighting the system and reducing thermal storage capacity to unacceptably low levels.
The present study addresses the impact of employing a layer of a porous medium over the inner hot pipe of an LHTES unit in the presence of nanoadditives. The porous layer can be placed with eccentricity toward the bottom to keep the advantage of natural convection flows at the top region of the TES unit while enhancing the conduction heat transfer at the bottom solid regions. The degree to which a metal foam layer affects the response (melting) time and the effect of nanoparticle additions on the heat transfer rate are investigated. The effectiveness of nanoparticles in the presence of the convective flows and the metal foam and the melting behavior of such a system will also be considered. A mathematical model will be developed to consider these effects. Figure 1a illustrates the structure of a TES unit utilized in this work. The inner pipe in the system is surrounded by a porous annulus with a permeability, κ, porosity, ε, and eccentricity, e. The velocity of a heat transfer fluid flowing in the inner pipe is high. In addition, the inner pipe employed in the system has a low thickness and a high thermal conductivity. Hence, it is reasonable to assume that the tube temperature is fixed at T h . The pores of the porous annulus and the void space between the porous annulus and the outer shell of the system are occupied with a nano-enhanced phase change material (NePCM) with homogenously dispersed nanoadditives of cooper or graphene oxide. The volume fraction of the nanoparticles in the base phase change material is σ. The outer shell of the system is well insulated, and the NePCM is initially at a supercooled temperature T c < T h . The geometrical characteristics of the cross-section of the unit (Figure 1b) are r s = 30 mm, r h = 2r s /5, and r p (e = 0) = 3r s /4. Table 1 lists the thermal properties of the PCM, the nanoadditives, and the metal foam. Table 1. Thermophysical characteristics of the materials used in this study [37,38].

Governing Equations
Liquid flow in the clear and porous domains can be solved by the controlling equations expressed below [39][40][41]: Continuity equation: ∂u ∂x Momentum equation in the x direction: Momentum equation in the y direction: where u and v are the velocity components for an NePCM, and t is the time; The symbols ρ and µ are the density and the dynamic viscosity, respectively; κ and ε are the permeability and porosity of the metal foam, respectively; s(T) is a source term which forces the velocities to zero in the solid domain; the subscripts of NeP and l indicate the NePCM and liquid state for the NePCM, respectively. In the above equations, spatial functions for the porosity and permeability ε j and κ j were defined as: where 1 and 2 in the above functions denote the porous and clear zones, respectively; κ is the permeability of the porous zone for a porous medium with a pore diameter, d p , and ligament diameter d l , was given by [42,43]: where Moreover, the momentum sink term, i.e., s(T), was defined as the following: where A mush is a large number O(10 6 ) and ϑ is a small number O(10 −3 ); η(T) is written as: where i describes the different states of the NePCM, C p is the heat capacity, and η is the liquid volume fraction; the subscripts of m and s denote the porous matrix and solid state of the NePCM, respectively. The thermal conductivity of the zones, λ m,j , was described as: where A large number of relationships which can be used to evaluate the effective thermal conductivity of the porous medium are presented in the literature [44]. In this study, the following relation, which takes into account the major characteristics of the porous matrix, was employed as following [28,45]: where

The Thermal and Physical Specifications of the PCM Enhanced with Nanoparticles
The density of the PCM enhanced with nanoparticles, ρ NeP , was expressed as: where σ is the concentration of the nanoadditives. The dynamic viscosity of the PCM enhanced with nanoparticles, µ NeP,1 , was written as [36]: The thermal-volume expansion coefficient of the PCM enhanced with nanoparticles, β NeP,1 , is expressed as: The effective thermal conductivity of the enhanced PCM, λ NeP,1 , was expressed as: The heat capacity of the PCM enhanced with nanoparticles, C p,NeP , was described as: the heat latent of the composite material, h sf,NeP , was written as:

Boundary Conditions
The controlling boundary conditions are as following: (a) at the interface of the porous and clear zones: (b) at the inner tube and shell: where T h = 316 K. The initial temperature was also considered as 304 K which is one degree below the fusion temperature.

Characteristic Parameters
The energy stored in the system with the nanoparticle-enhanced PCM, ES, was defined by: Finally, the melting volume fraction, MVF, was evaluated as:

Solution Approach
The system is divided into three regions: the first region where the PCM is completely solid, the second region where it is completely molten, and the third region where solid and liquid states of the material coexist (the mushy zone). In order to distinguish these three zones, the fusion temperature and the melting temperature window, δT, are used. Melting and solidification can be described by a movement of the interfaces between these three zones and a change in the phase fraction of liquid within the mushy zone. The melting process can be modeled by using source terms for momentum, s(T), and energy, via the rate of latent heat generation and is proportional to the rate of change of fraction of liquid, ∂η(T)/∂t, which are described in the previous section. In addition, the velocity components are controlled by the source term, s(T). This source term acts as a virtual porous material with variable permeability, as introduced in Equation (5). The magnitude of s(T) approaches zero in the liquid domain (no force on the fluid) and has a very large number (impermeable medium) in the solid domain. Since the velocity of liquid reaches zero in the solid domain through the mushy zone, a high-velocity gradient is seen inside the mushy zone. The energy source term, ∂η(T)/∂t, acts as a heat sink within the mushy zone, controls the latent heat abortion and results in a high-temperature gradient. To calculate the high gradients in the mushy zone, a very fine mesh is vital. However, using such a fine mesh throughout the domain would be prohibitive in terms of computational expense. Therefore, a mesh adaptation technique was used.
A Galerkin finite element method technique was used to solve the nonlinear differential equations in this study. The phase change equations were implemented by user defined codes. According to this method, the equations and the corresponding boundary and initial conditions were transformed to a weak form and solved by the Galerkin finite element method. The residual equations based on shape functions were introduced and solved iteratively. The numerical approach could be found in [46,47], and it was not repeated here for the sake of brevity.
A dummy phase field variable, η 0 , was introduced that acted as a criterion for grid adaption, i.e., only when η 0 = 1, the mesh adaptation will occur. η 0 was defined as: In order to control the time steps, a free step backward differentiation formula with an automatic time step range of 1-2 was used [48]. The Newton method employing PARallel DIrect SOlver (PARDISO) solver was adopted [49][50][51]. A residual error of~10 −6 and a Newtonian damping factor of 0.8 were considered to solve the residual equations iteratively.

Mesh Sensitivity Analysis
Five different meshes were considered for testing the mesh sensitivity according to Table 2. As seen in Figure 2, the MVF and the melting interface after 500 s for the PCM free of nanoparticles (ε = 0.80, e = 0.22r s , and σ = 0.00) were used. There is a good agreement in all cases, especially cases III, IV, and V. Case III is chosen for further studies, as it is the least computationally expensive. There are 5104 elements in the domain, of which 335 lie at the boundaries. The mesh deformation and the adaptive refinement at t = 0 and t = 500 are demonstrated in Figures 3 and 4, respectively. A desirable transient can be seen. Moreover, the results of the mesh sensitivity analysis showed a negligible error by using case III. Thus, the numerically uncertainty is negligible and will be not shown in the reported results.

Validation with the Literature
Previous experimental and numerical investigations have been used to validate the results and check the correctness of the numerical simulation of the present study. An experimental neutron flux study of the melting process and the evolution of the interface between the solid and liquid phases under a constant heat flux on the vertical walls of a square enclosure and other insulated walls gave an interface that has the same shape as is expected from the output of the simulation proposed in this study ( Figure 5; [52]). A second study into the phase change of a paraffin wax in a copper foam with a porosity of 0.95 was compared to the output of the proposed simulation at three different time steps, and the shape of the interfaces was found to be similar in all three cases (Figure 6; [45]). A third study was also used, in which a pure (free of nanoparticles) PCM was melted in a square enclosure with the two side walls each held at a separate, constant temperature to establish a constant thermal gradient across the material. The results of that study were in agreement with the results of the model proposed here (Figure 7; [53]) at two non-dimensional times of Fo = 2.05 and 3.48. Finally, the calculated heat transfer by natural convection was compared to experimental results of the zone between two horizontal cylinders. The isotherm contours of two studies were compared at Pr = 6.19 and Ra = 2.33 × 10 5 and show reasonable agreement (Figure 8; [54]).

Results and Discussion
The model was used to investigate the effects of eccentricity e within the range of 0-0.22 r s , porosity ε within the range of 0.8-1.0 (i.e., the metal foam occupies a volume between 20% and 0% of the available space, such that a porosity of 1 implies that no foam is present), the type of nanoparticles added, and the volume fraction σ of each type of nanoparticles ranging from 0.00 to 0.08.
When there is no foam (ε = 1), a layer of the melted PCM remains around the inner hot cylinder, while for the other values of ε, the PCM melts near the hot wall at first but the molten region then expands and covers the whole annulus ( Figure 9). This implies that the presence of the copper foam contributes significantly to the melting of the PCM. Heat is transmitted by conduction within the foam, and this causes the PCM to melt. Once the PCM melts, convection in the melt is also able to transfer heat: hot liquid rises while cool liquid falls. Recirculation zones appear in the melt and surround the inner cylinder, which contributes further to melting.
In the case of a concentric annulus (no eccentricity, e = 0), the recirculation zone starts symmetric around the inner cylinder ( Figure 10). As time progresses, convective effects become more pronounced and more portions of the PCM melt above the inner pipe and some remain solid below the inner tube. When eccentricity is increased, the porous zone around the inner pipe is shifted downwards. As the presence of the foam promotes melting as indicated earlier, all the PCM portions below the inner tube eventually melt. In the upper part of the cavity, the dominant convective effects also lead to full melting of the PCM in that region.
For a foam-free system (ε = 1), the melt volume fraction increases slowly with time. For the other values of ε, the melt volume fraction increases much faster. This is because a lower value of ε implies a greater presence of the solid foam and, consequently, more heat transmission throughout the PCM, which ultimately leads to melting (Figure 11). It is also shown that for a concentric system (e = 0), the time required to achieve complete melting is higher, compared to that in an eccentric system. This is related to the observations shown in Figure 10. When there is no eccentricity, convection enhances melting in the upper region where hot liquid is ascending, while at the bottom of the annulus, melting is slow as the copper foam does not fully cover that region. When eccentricity is greater, the PCM melts faster in the bottom region, while melting also continues in the upper region due to free convection ( Figure 12). The full melting time is approximately halved in a system with an eccentricity equivalent to 0.1r s , compared to in a concentric system.    For a greater fraction of nanoparticles (increased σ), the melting front is more developed and closer to the outer cylinder, indicating an increase in the molten PCM. This is true for both copper and graphite oxide nanoparticles. The nanoparticles dispersed in the PCM improve heat transfer. The type of nanoparticles has a limited impact on the melting front, as the location of the front remains almost the same for both graphite oxide and copper nanoparticles, as does the volume fraction of the PCM ( Figure 13). As seen, using 8% of nanoparticles could reduce the melting time by 13%, caculcated by 100 × (1350−1175)/1350 for both types of nanoparticles. Initially, the volume fraction of the molten PCM increases sharply, as the conduction in the copper foam is dominant at short times. The volume fraction of the molten material then increases gradually when convection dominates, until melting is complete. In all cases, the time to complete melting is lower, when nanoparticles are present. The melting time is minimized, when the amount of nanoparticles is greatest (σ = 8%). Total melting occurs faster when copper nanoparticles are used compared to that when graphite oxide particles are used.
The latent storage capacity and charging power are greater than the sensible capacity and power (Figure 14), as most of the involvement of the PCM in the heat transfer is consumed by the melting process. It is shown that the latent and sensible storage heat capacities are almost the same for all the values of eccentricity since the phase change mass is constant. Since the stored heat is constant, the charging power is only sensitive to the melting rate. Melting is accelerated due to the presence of the porous foam around the inner pipe and convection in the upper region of the cavity, and e = 0.1r s corresponds to the optimal placement of tube under those conditions, compared to zero eccentricity. Conversely, when the inner pipe and the surrounding porous zone are concentric, the melting time is greatest, and consequently, the average charging power is lowest.
Reducing porosity (raising ε) leads to a slight reduction in the sensible storage capacity and an increase in the latent capacity ( Figure 15). The average charging power decreases sharply to almost zero, when there is no foam (ε = 1). In this case, the total melting time is very long, as no heat is transmitted by the metal foam, which is consistent with the data in Figures 9 and 11. A moderate decrease in the storage capacity is observed, when the nanoparticle content σ is increased for both copper and graphite oxide ( Figure 16). Charging power is higher when copper nanoparticles are used, compared to when graphite oxide nanoparticles are used. This is consistent with the observation that the total melting time is lower in the case of copper nanoparticles and when a lower volume fraction of nanoparticles is used ( Figure 13).

Conclusions
A shell-and-tube LHTES system was modeled, and its performance was simulated. The inner tube was supported by a layer of copper metal foam. Graphite oxide and copper nanoparticles were tested as additives to accelerate heat transfer in a PCM. The impact of convective heat transfer in the molten region was taken into account. The outcomes of the present numerical analysis can be summed up as follows:

•
The presence of a metal foam around the inner pipe contributes significantly to heat transfer and to the melting of the PCM. Reducing the porosity of the foam increases the amount of heat transmitted through the foam and enhances the melting of the PCM.

•
The eccentricity of the system affects heat transfer by changing the location of the foam around the inner pipe. The melting of the PCM is limited by convection in the upper part of the cavity and by conduction through the metal foam in the lower region of the cavity. The melting of the PCM is slowest, when the inner pipe and the porous zone are concentric; the total melting time is halved for the eccentricity e = 0.1r s . • Using a higher volume fraction, σ, of nanoparticles enhances melting. Copper nanoparticles are slightly more effective in transferring heat than graphite oxide nanoparticles. For both types of nanoparticles, the time required to melt the PCM completely is 13% shorter for an 8% volume fraction of nanoparticles (σ = 0.08).

•
The three parameters, the porosity ε, the eccentricity e, and the volume fraction of nanoparticles, σ, have limited effects on the total storage capacity. However, the average charging power and the total melting time are affected more strongly by these parameters. Overall, the charging power decreases for longer total melting time, mainly for smaller values of eccentricity e and for lower volume fractions of foam (higher values of ε). The charging power is highest, when eccentricity takes the value e of 0.1r s .

Conflicts of Interest:
The authors declare no conflict of interest.