Thermal Charging Optimization of a Wavy-Shaped Nano-Enhanced Thermal Storage Unit

A wavy shape was used to enhance the thermal heat transfer in a shell-tube latent heat thermal energy storage (LHTES) unit. The thermal storage unit was filled with CuO–coconut oil nano-enhanced phase change material (NePCM). The enthalpy-porosity approach was employed to model the phase change heat transfer in the presence of natural convection effects in the molten NePCM. The finite element method was applied to integrate the governing equations for fluid motion and phase change heat transfer. The impact of wave amplitude and wave number of the heated tube, as well as the volume concertation of nanoparticles on the full-charging time of the LHTES unit, was addressed. The Taguchi optimization method was used to find an optimum design of the LHTES unit. The results showed that an increase in the volume fraction of nanoparticles reduces the charging time. Moreover, the waviness of the tube resists the natural convection flow circulation in the phase change domain and could increase the charging time.


Introduction
The major threat in energy society is the expanding gap between the supply of energy and global demand, which makes it challenging for engineers to find unique solutions to fulfill the needs of this society. The problems on energy conservation and storage were under the scientific community's concern since figuring out the idea of utilizing the industrial hot waste streams and the abundant-solar thermal energy [1]. The latent heat is an alternative to store thermal energy in the form of latent heat of phase change through a melting or freezing process [2]. The classifications of phase change materials (PCMs) include eutectics, organic, and inorganic materials. Most studies focus on organic materials due to their large heat storage capacity, different phase-change temperatures, and low cost [3,4]. A PCM for thermal energy storage was prepared by Osterman et al. [5] for heating and cooling purposes to save energy. Their study showed annual save reaches 142 kWh in the energy consumption of an office. Substitutes for PCM have been developed recently by the composite metal foam/PCM to enhance the heat transfer inside PCM [6,7]. A study by Mesalhy et al. [8] showed that the porous foam matrix within PCM has a significant impact on the rate of heat transfer and charging time. Sardari et al. [2] made a numerical simulation to evaluate the discharging mechanism in a phase change material to air heat exchanger, which aimed to have space heating using PCM and composite of copper foam. Significant advantages were shown in the result. By using the composite PCM, a 73% higher heat retrieval rate was achieved, and the solidification time was reduced by 43%. The low thermal conductivity is considered a well-known property of organic PCMs. This property leads to a limitation in absorbing and releasing heat. Therefore, nanoparticles provide good physical and chemical properties to the PCMs compared to their bulk form [9,10]. Achieving higher thermal conductivity is related to the reduction in nanoparticle size [11]. A new sort of nano-enhanced phase change material (NePCM) was developed by Wu et al. [12] by suspending a small amount of C/Cu, Cu, and Al nanoparticles paraffin to improve the heat transfer rate of paraffin. The study found that Cu nanoparticles have the best heat transfer performance.
Chen et al. [13] added a small amount of CuO nanoparticles to paraffin to improve thermal storage and solar thermal conversion. The nanoparticles impressively improved the light abortion and solar energy conversion. The nanoparticles could improve the steady temperature of the unit up to 2.3 times compared to the pure paraffin.
Developing additives is an interesting task that has been advanced by the synthesis of novel nanoadditives. Successful utilization of NePCM in heat transfer applications depends on some important terms, which can be summarized as follows [14]: (1) on the level of required flexibility, thermal stability, strength, and corrosion resistance. (2) Protection shield for the PCM to avoid any harmful interaction with the environment. (3) Make a sufficient heat transfer surface. (4) Provide easy handling and structural stability.
One of the effective and economical ways for heat transfer enhancement is improving the heat transfer area by using fins and extended surfaces. However, fins could obstruct the free convection circulation in molten regions [15]. Hence, a smooth development of convection by employing fins can only be achieved by optimized fin size [16]. A combination of using nanoparticles and an optimum structure of fin could avoid the disadvantage of employing nanoparticles and fins in a conical shaped latent heat thermal energy storage (LHTES) system. Singh et al. [15] used a combination of active (fins and nanoparticles) heat transfer enhancement techniques to make an innovative LHETS, where the volume of the fins was fixed (0.15) for all cases. The study found a 57% reduction in charging time of LHTES with a conical shape in the presence of 5% graphene nano-plates. On the other side, no substantial improvement was observed in thermal performance of finned LHTES with 1% GNP (graphene nanoplatelets) because of the restriction of natural convection vortices and raised viscosity. Different studies found that the geometry of PCMbased heat storage units influences heat transfer performance such as the number and the shape of heat transfer fluid (HTF) inner tubes [17]. Mehta et al. [18] utilized different tubes in an LHTES, i.e., multiple-tubes, the eccentricity of heat transfer fluid (HTF) tube, and rotation of the HTF tube. The study concluded that the maximum reduction in the melting time was achieved by rotating latent heat storage (LSH) systems with a 47.75% reduction. Wavy walls are considered a promising technique on heat transfer enhancement instead of straight walls by mixing up the flow in both sides of the wavy wall and increasing the surface area of heat transfer [19][20][21][22]. Shahsavar et al. [23] investigated the solidification and melting performance in a wavy LHTES with composite PCM and copper foam. They showed that a wavy enclosure could have promising advantages in the system, where strong influences of metal foams suppressed the convection heat transfer.
The literature studies show the impact of wavy walls on the convective and conductive heat transfer flows. Thus, the present study aims to address the impact of using wavy tubes in shell and tube heat thermal energy storage for the first time. The current research aims to investigate the advantage of using wavy tubes in clear flows (with no porous metal foam) where the convective heat transfer effects are important. The impact of using nanoadditives on the charging time of the energy storage system will also be investigated.

Mathematical Model
A multi-tube shell-and-tube LHTES unit, in which the NePCM is arranged in the shell side with an initial super cold temperature of T initial , and the hot fluid flows through the wavy tubes, as shown in Figure 1a. The tubes are placed next to each other. Here, one of the tubes (see Figure 1b) is selected for the investigation of melting heat transfer. The HTF with a high temperature is injected into the wavy tube with a temperature of T in , the heat is transferred between the wavy wall and the NePCM, with thermal energy stored in the NePCM. A 2D axis symmetric view of the model is shown in Figure 1c.
The literature studies show the impact of wavy walls on the convective and conductive heat transfer flows. Thus, the present study aims to address the impact of using wavy tubes in shell and tube heat thermal energy storage for the first time. The current research aims to investigate the advantage of using wavy tubes in clear flows (with no porous metal foam) where the convective heat transfer effects are important. The impact of using nanoadditives on the charging time of the energy storage system will also be investigated.

Mathematical Model
A multi-tube shell-and-tube LHTES unit, in which the NePCM is arranged in the shell side with an initial super cold temperature of Tinitial, and the hot fluid flows through the wavy tubes, as shown in Figure 1a. The tubes are placed next to each other. Here, one of the tubes (see Figure 1b) is selected for the investigation of melting heat transfer. The HTF with a high temperature is injected into the wavy tube with a temperature of Tin, the heat is transferred between the wavy wall and the NePCM, with thermal energy stored in the NePCM. A 2D axis symmetric view of the model is shown in Figure 1c. The PCM and the nanoadditives on the NePCM are the coconut oil and CuO nanoparticles, with the thermophysical characteristics listed in Table 1. The length of the system, i.e., L, is 100 mm, the outer radius of the computational domain, i.e., Ro, is 50mm, and the radius of the simple tube, i.e., Ri, is 9 mm. Assuming the tube is sufficiently thin and high-thermally conductive, the temperature gradient and energy stored can be ignored in the tube. The amplitude and period of the wavy wall are A and N, respectively. The properties of the HTF flowing in the tube, i.e., water, are tabulated in Table 1. 103 ± 1% --- The PCM and the nanoadditives on the NePCM are the coconut oil and CuO nanoparticles, with the thermophysical characteristics listed in Table 1. The length of the system, i.e., L, is 100 mm, the outer radius of the computational domain, i.e., R o , is 50 mm, and the radius of the simple tube, i.e., R i , is 9 mm. Assuming the tube is sufficiently thin and high-thermally conductive, the temperature gradient and energy stored can be ignored in the tube. The amplitude and period of the wavy wall are A and N, respectively. The properties of the HTF flowing in the tube, i.e., water, are tabulated in Table 1. Table 1. Properties of LHTES materials [24,25].

Convective Phase Change Heat Transfer in NePCM
The NePCM melts over time, and the natural convection of the molten NePCM occurs in the liquid phase. The transient phase change process with natural convection is described using the volume-averaged approach: (I) Mass conservation 1 r ∂(ru r ) ∂r where r and z are radial and axial directions while u r and u z are radial and axial velocity components.
(II) Momentum equations where g, µ, ρ, and β are the gravity acceleration, dynamic viscosity, density, and thermal expansion coefficient, respectively. The filed variables of p and T are the pressure and temperature, respectively while t is time. The subscripts l, NP, and me refer to liquid, NePCM, and fusion temperature, respectively. NePCM. Terms containing A mush are applied to damp the velocity vector in the solid phase. The numerical constant A mush and λ are relatively large (5 × 10 5 kg/(m 3 s)) and small (10 −3 ), respectively. The function γ(T) of the above equations is described as the following: where T is temperature and ∆T me is the melting range.
(III) Energy conservation in which k, h f , C p , and ω na are the thermal conductivity, specific heat capacity, latent heat, and volume concentration of nanoparticles, respectively. The subscripts PCM and na show the pure PCM and nanoparticles. In the above equations, the property of NePCM is the weighted average of its liquid and solid states, depending on the melt fraction γ(T), which are evaluated as: where subscript s indicates the solid state.

The Properties Model
The NePCM is a mixture of PCM and nanoparticles. Thus, following the conservation of mass, the density can be evaluated as: in which ω na denotes the volumetric concertation of CuO nanoparticles. The Brinkman model is applied to evaluate the liquid NePCM's dynamic viscosity as: The volumetric expansion of NePCM is important for the evaluation of buoyancy forces. This coefficient for the NePCM mixture is evaluated as: The Maxwell relation is employed to estimate the NePCM's thermal conductivity as [26]: Equation (9) was discussed in [26] and showed good agreement with experimental measurement for 15-25 nm spherical nanoparticles. Finally, the NePCM's heat capacity was computed as:

Convective Heat Transfer in HTF
Assuming the HTF injected in the tubes is laminar, transient, and incompressible, the convective heat transfer can be described by the below equations: (III) Energy conservation where the subscript HTF refers to heat transfer fluid.

Controlling Boundary and Initial Conditions
The temperature and heat continuity are adopted for the thermal boundary condition at the surface of the wavy tube as [27,28]: The HTF is injected with the temperature of T in and v in .
where T in is the inlet temperature for HTF. At the outlet of the HTF tube, the developed convective flow is considered to be established.
Due to the symmetrical distribution of the tube bundle, the boundary conditions of the computational domain are: At the upper and lower surfaces of the tube bundle, the boundary conditions can be defined as the following: The initial condition for the NePCM domain can be expressed as the following: where in the above equation T in = 323 K, and T initial = 293 K.

Heat Transfer Characteristics
The cumulative heat transferred (ES(t)) from the HTF to the NePCM is exactly the stored energy (ES) in the unit. ES(t) can be calculated by integrating the instantaneous heat transfer rate from the wavy wall: In Equation (21), the first term shows the sensible heat stored in the NePCM, and the second term takes into account the latent heat stored in the NePCM. Moreover, the rate of the absorbed heat by NePCM, ∆q, can be computed as: Moreover, the average temperature of the interface surface can be evaluated as: The average convective heat transfer coefficient (h HTF ) at the water side can be computed by using Newton's cooling law as: ∆q = 2πR i × L × h HTF × T in − T s where the inlet temperature and the enclosure high were considered as the reference temperature and the characteristic length, respectively. Thus, h HTF is computed as: Using the same approach, the heat transfer coefficient (h NP ) at the NePCM side can also be computed as following. By neglecting small unsteady temperature fluctuations at the water side, the transported energy from the water to the NePCM is equal to stored energy in the NePCM. Thus, h PCM can be introduced as: Here, L was taken as the characteristic length of the wavy tube instead of the arc length s. This is because the variation of amplitude and the wave number of the wavy pipe could change the arc length, and hence, the using arc length s would lead to the convection heat transfer directly dependent on the multiplication of heat transfer coefficient and the actual surface of heat transfer. However, since L is constant, the computed values of h can be directly compared and judged. Moreover, the thermal resistance can be computed as: Furthermore, the pressure loss of the HTF passing the channel can be expressed as the following: the total Melting Volume Fraction (MVF) is evaluated as: In the practical application of the LHTES unit, the performance of the unit can be evaluated by charging power. This parameter shows the capacity of the PCM to store energy and depends on the amount of the energy stored at 100% melting volume fraction:

CP =
ES when the melting process complete complete melting time ,

Numerical Approach and Mesh Study and Verifications
Since the configurations of a single tube along with the surrounding NePCM is axisymmetric, a two-dimensional (2D) configuration can be considered to save computational time, as shown in Figure 1c. As previously discussed, in detail, to model the melting front, the enthalpy-porosity method is employed. The time-dependent finite element method involving the Galerkin approach is employed to solve the temperature, pressure, and velocity fields. Convergence is obtained for each simulation when scaled residuals are less than 10 −3 for mass, momentum, and energy conservation equations. The influence of mesh size on the accuracy of the results is examined through a mesh sensitivity analysis for MVF. A structured mesh due to efficiency and simplicity is used to discretize the domains. Naturally, an unstructured grid requires more memory than a structured grid with an equal number of elements. As tabulated in Table 2, the liquid fractions of NePCM for five grids with different elements count are compared. Through the relative errors, it is found that a mesh with a size of 30 × 150 elements can be adequate for the purpose of optimizing the mesh elements. The results of the model discussed above are compared with the numerical and experimental results reported in the literature [29][30][31][32][33]. Kumar et al. [29] experimentally and numerically studied the fusion process of lead contained in a cuboid exposed a constant heat flux at the right sidewall. The dimensions of the cuboid were 50 mm in height, 50 mm in width, and 60 mm in depth. Figure 2 illustrates the progress of the melting front in the experimental and numerical study conducted by Kumar et al. [29] and the present numerical model solved by FEM. Gau and Viskanta [30] experimentally analyzed the buoyancy-driven flow in the melted pure PCM and its impact on the melting front motion. In this experiment, the melting process was followed in a rectangular test cell with a heat source located on the left sidewall and a heatsink on the opposite sidewall. The dimensions of the test cell used were 63.5 mm in width, 88.9 mm in height, and 38.1 mm in depth. Figure 3 depicts a set of progressive solid-liquid interfaces of the numerical studies compared with the experimental interfaces observed by Gau and Viskanta. There are good agreements between the results of the melting front location in the current study and the results of other works [29][30][31][32][33]. Hence, the utilized code is reliable for parametric investigation. Although it would not provide the exact values of the physical quantities, it can provide the effect of each parameter on the system. The results of the model discussed above are compared with the numerical and experimental results reported in the literature [29][30][31][32][33]. Kumar et al. [29] experimentally and numerically studied the fusion process of lead contained in a cuboid exposed a constant heat flux at the right sidewall. The dimensions of the cuboid were 50 mm in height, 50 mm in width, and 60 mm in depth. Figure 2 illustrates the progress of the melting front in the experimental and numerical study conducted by Kumar et al. [29] and the present numerical model solved by FEM. Gau and Viskanta [30] experimentally analyzed the buoyancy-driven flow in the melted pure PCM and its impact on the melting front motion. In this experiment, the melting process was followed in a rectangular test cell with a heat source located on the left sidewall and a heatsink on the opposite sidewall. The dimensions of the test cell used were 63.5 mm in width, 88.9 mm in height, and 38.1 mm in depth. Figure 3 depicts a set of progressive solid-liquid interfaces of the numerical studies compared with the experimental interfaces observed by Gau and Viskanta. There are good agreements between the results of the melting front location in the current study and the results of other works [29][30][31][32][33]. Hence, the utilized code is reliable for parametric investigation. Although it would not provide the exact values of the physical quantities, it can provide the effect of each parameter on the system.

Taguchi Analysis
Taguchi's technique is a powerful statistical method utilized to optimize designs, in which the results are influenced by different parameters or control factors. Each of these control factors is assigned a number of numerical values at different levels. Taguchi's method is then applied in order to find a minimal number of test cases required to determine the optimal values of the control factors. This is more efficient than performing all the possible tests, which presents a high cost in time and resources and can be practically very difficult to implement.

Taguchi Analysis
Taguchi's technique is a powerful statistical method utilized to optimize designs, in which the results are influenced by different parameters or control factors. Each of these control factors is assigned a number of numerical values at different levels. Taguchi's method is then applied in order to find a minimal number of test cases required to determine the optimal values of the control factors. This is more efficient than performing all the possible tests, which presents a high cost in time and resources and can be practically very difficult to implement.
In the present study, the key parameters of the problem are the volume fraction of the nanoadditives (0 ≤ ω na ≤ 0.05), the amplitude of the wavy wall (0 ≤ A ≤ 3mm), the period of the wavy wall (1 ≤ N ≤ 5). These parameters represent the control factors, and each is assigned four levels. Table 3 summarizes the control factors of the problem and their correspondent levels. If all the possible combinations of the control factors and their levels were to be tested, 4 3 tests would be required. The L16 orthogonal array is then used, and Taguchi's algorithm is applied in order to reduce to16 the number of required tests. The range of the control factors and their levels corresponding to these 16 tests are listed in Table 4. In this study, the computation continued until the unit reaches full melting, i.e., MVF = 1. As previously mentioned, the aim of Taguchi's method is to determine the values of the control factors that would lead to an optimal result. In this regard, based on the desired result, one of the following three approaches can be used: the larger-the better, the nominal-the better, and the smaller-the better. In the present study, the time required for complete PCM melting is considered as the optimization result. As the objective is to reduce the full melting time, then the smaller-the better approach is adopted. The signal to noise ratio, SNR, is analyzed for each test. SNR represents the advantage of wanted or desired results to the undesired results. Moreover, a linear regression model is developed in order to relate the full melting time to the control factors: After determining the combination of parameters that lead to the lowest full melting time among the tests, the optimal values of the control factors can be defined. These values are presented in Table 5: ω na = 0.06 (level 4), A = 0 (level 1), and N = 4 (level 4) with a full melting time of 22,709 s. Since A = 0, the value of N is meaningless, and indeed, the Taguchi method shows that a plane tube with no wave could lead to better performance compared to a wavy tube. The next step is to define which control factors have the most influence and the results. To do so, the SNR of the tests are analyzed. The mean values of the SNR for each test are plotted in Figure 4, while Table 6 list the values of ∆ indicating the difference between the maximum and minimum of the mean SNR in each case. It is evident from Figure 4 that the optimal values previously obtained lead to the highest mean value of SNR among the tested parameters. ∆ is an indicator of the degree of influence of each control factor, i.e., the higher the value of ∆, the greater the influence of the corresponding parameter on the target-parameter. Thus, based on Table 6, the various control factors can be classified in the order of their influence: A > ω na > N.   To further confirm the design optimization with the previously determined values, nine more tests are conducted with combinations of the design parameters around the optimal values. In each group test, two of the three control factors are kept constant at their optimal value, while the third is varied among three different values. The results show that the full melting time in all the test cases is greater than 22,709 s, which corresponds to the optimal values. This confirms that the design is optimized for the values of the control factors previously found: ωna = 0.06, A = 0, and N = 4.  To further confirm the design optimization with the previously determined values, nine more tests are conducted with combinations of the design parameters around the optimal values. In each group test, two of the three control factors are kept constant at their optimal value, while the third is varied among three different values. The results show that the full melting time in all the test cases is greater than 22,709 s, which corresponds to the optimal values. This confirms that the design is optimized for the values of the control factors previously found: ω na = 0.06, A = 0, and N = 4.

Parametric Investigation
As the optimal values are determined, the next step is to analyze the flow and thermal patterns in the domain in order to provide physical insights explaining the aforementioned outcomes. Table 7 shows nine selected designs around the optimum case. Cases 1-3, 4-6, and 7-9 aim to investigate the impact of concentration of nanoparticles, wave amplitude, and wave number on the melting process, respectively. Figures 5 and 6 depict the isothermal contours and the streamlines in the cavity for two values of the nanoparticle volume fraction, ω na = 0.06 corresponding to the optimal case and ω na = 0 corresponding to pure PCM without nanoparticles. The presence of streamlines indicates that the PCM is in the liquid phase; i.e., it has melted in that region. Initially, the PCM starts to melt near the heated tube at the left. In that region, the conductive effects are important. As time goes, the hot melted PCM circulates upwards due to the density difference, and free convection takes place.
A clockwise vortex is created in the domain. The hot PCM impinging in the top region transfers its heat to the melting front so that the convective effects, and consequently PCM melting, are very weak in the bottom part of the cavity compared to its upper part. It can be seen that PCM melting occurs faster in the case ω na = 0.06, which is due to conduction heat transfer enhancement at the early stage of melting. It should be noted that the presence of nanoparticles increases the dynamic viscosity of liquid PCM and reduces the strength of natural convection flows at the final stages of melting. This observation is further confirmed in Figure 7, which illustrates the evaluation of the melted volume fraction MVF and the stored energy ES as functions of time for different values of ω na . In the charging process, both the MVF and the ES increase for higher values of ω na . Moreover, full melting (MVF = 1) is achieved faster when ω na is increased. This is due, as previously explained, to the enhanced thermal conductivity when nanoparticles are added to the PCM, and the growth of dynamic viscosity. The stored energy at 750 min is slightly larger for pure PCM compared to NePCM, which is mainly due to the higher overall latent heat of pure PCM. The presence of nanoparticles reduces PCM's overall latent heat since they cannot store latent heat. They also slightly increase the heat capacity of the PCM. However, the increase of heat capacity contributing to the sensible heat is much lower than the change in the latent heat. Thus, by using nanoparticles, the ultimate value of stored energy can be reduced slightly. Figures 5 and 6 depict the isothermal contours and the streamlines in the cavity for two values of the nanoparticle volume fraction, ωna = 0.06 corresponding to the optimal case and ωna = 0 corresponding to pure PCM without nanoparticles. The presence of streamlines indicates that the PCM is in the liquid phase; i.e., it has melted in that region. Initially, the PCM starts to melt near the heated tube at the left. In that region, the conductive effects are important. As time goes, the hot melted PCM circulates upwards due to the density difference, and free convection takes place.  A clockwise vortex is created in the domain. The hot PCM impinging in the top region transfers its heat to the melting front so that the convective effects, and consequently PCM melting, are very weak in the bottom part of the cavity compared to its upper part. It can be seen that PCM melting occurs faster in the case ωna = 0.06, which is due to conduction heat transfer enhancement at the early stage of melting. It should be noted that the presence of nanoparticles increases the dynamic viscosity of liquid PCM and reduces the strength of natural convection flows at the final stages of melting. This much lower than the change in the latent heat. Thus, by using nanoparticles, the ultimate value of stored energy can be reduced slightly. Figure 7c displays the thermal resistance (RNP) at the NePCM side for various nanoparticle concentrations. Generally, RNP increases as the melting advances. By the advancement of melting, the distance of the cold surface (melting front) increases from the heated surface (wavy surface). The increase of distance increases the thermal resistance. This figure shows that at the beginning of melting, where the heat transfer is conduction dominant, the increase of volume fraction of nanoparticles decreases the thermal resistance, RNP. About 300 min, there is a turning point that the thermal resistance grows as the volume fraction of nanoparticles rises. Figure 7a shows about 80% melting at 300 min. Thus, at 300 min, a significant amount of liquid PCM exists in the enclosure, and the natural convection heat transfer is significant (see Figure 6). Therefore, nanoparticle presence increases the dynamic viscosity and reduces the strength of natural convection leading to the growth of thermal resistance.
The effect of the hot tube amplitude on the temperature distribution and on the flowpatterns is illustrated in Figures 8 and 9. Two amplitudes are considered: A = 0 (flat tube), which corresponds to the optimal case, and A = 0.06. It is clear that the streamlines mimic the wavy tube profile in the left portion, where heat is mainly transmitted by conduction. This behavior remains when time passes. However, it is more apparent in the bottom part of the cavity, as convective effects dominate the PCM behavior in the upper part. For Cases 4-6, the steady-state pressure drop was computed as 2.01 × 10 −4 Pa (A = 1 mm), 6.63 × 10 −3 Pa (A = 2 mm), and 6.53 10 −2 Pa (A = 3 mm). The pressure drop of the optimum case is 3.55 × 10 −5 Pa (A = 0 mm). As seen, there is an exponential pressure drop by the increase of undulation amplitude.  Figure 7c displays the thermal resistance (R NP ) at the NePCM side for various nanoparticle concentrations. Generally, R NP increases as the melting advances. By the advancement of melting, the distance of the cold surface (melting front) increases from the heated surface (wavy surface). The increase of distance increases the thermal resistance. This figure shows that at the beginning of melting, where the heat transfer is conduction dominant, the increase of volume fraction of nanoparticles decreases the thermal resistance, R NP . About 300 min, there is a turning point that the thermal resistance grows as the volume fraction of nanoparticles rises. Figure 7a shows about 80% melting at 300 min. Thus, at 300 min, a significant amount of liquid PCM exists in the enclosure, and the natural convection heat transfer is significant (see Figure 6). Therefore, nanoparticle presence increases the dynamic viscosity and reduces the strength of natural convection leading to the growth of thermal resistance.
The effect of the hot tube amplitude on the temperature distribution and on the flowpatterns is illustrated in Figures 8 and 9. Two amplitudes are considered: A = 0 (flat tube), which corresponds to the optimal case, and A = 0.06. It is clear that the streamlines mimic the wavy tube profile in the left portion, where heat is mainly transmitted by conduction. This behavior remains when time passes. However, it is more apparent in the bottom part of the cavity, as convective effects dominate the PCM behavior in the upper part. For Cases 4-6, the steady-state pressure drop was computed as 2.01 × 10 −4 Pa (A = 1 mm), 6.63 × 10 −3 Pa (A = 2 mm), and 6.53 10 −2 Pa (A = 3 mm). The pressure drop of the optimum case is 3.55 × 10 −5 Pa (A = 0 mm). As seen, there is an exponential pressure drop by the increase of undulation amplitude.
For a fixed N and A, the values of R HTF were changed minimally during the phase change. This is since the HTF fluid-flow is independent of the melting process. The slight changes are from the local variation of the average temperature of the tube during the phase change. The average values of R HTF were found as 0.63 K/W (A = 1 mm), 0.73 K/W (A = 2 mm), 0.65 K/W (A = 3 mm). The R HTF resistance for the optimum design (straight tube) was 0.5 K/W, which is smaller than all wavy cases.  For a fixed N and A, the values of RHTF were changed minimally during the phase change. This is since the HTF fluid-flow is independent of the melting process. The slight changes are from the local variation of the average temperature of the tube during the phase change. The average values of RHTF were found as 0.63 K/W (A = 1 mm), 0.73 K/W (A = 2 mm), 0.65 K/W (A = 3 mm). The RHTF resistance for the optimum design (straight tube) was 0.5 K/W, which is smaller than all wavy cases. The waviness of the tube increases the distance between the heated line and the PCM and slows down its melting compared to a flat tube. Moreover, the waves act as obstacles that trap the convective flows in the molten region of the enclosure. Thus, the increase of surface contact is not as effective as the suppression of natural convection flow. This appears in the size of the thermal layer near the tube, which is greater for A = 0. In the upper part of the cavity, this behavior is less effective as the convective effects are dominant. On the other hand, PCM melts faster in the bottom part in the case of a flat-tube. This results in an overall faster PCM melting and a greater volume of liquid PCM for A = 0. For these reasons, it is clear that reducing A shows a positive effect on both MVF and ES as indicated in Figure 10. for A = 0. For these reasons, it is clear that reducing A shows a positive effect on both MVF and ES as indicated in Figure 10. The last control factor to be assessed is the number of undulations N of the wavy tube. The thermal and flow patterns of the PCM in the enclosure are shown in Figures 11  and 12 for two values of N, N = 4 (the optimal value) and N = 1, and for the highest wave amplitude of 3 mm. The flow streamlines mimic the wall profile near the hot tube, as the number of oscillations of the streamlines is equal to that of the tube. The volume of melted PCM, as well as the temperature in different zones of the cavity, are larger for N = 4 compared to N = 1, which can be attributed to the increase of surface contact between the PCM and the hot tube when the same amplitude is considered. This results in an increase of the MVF and the ES when N is raised, as observed in Figure 13. However, the presence of undulation entraps the fluid flow at the tube side and reduces the heat transfer. Thus, on the one hand, the presence of N contributes to the increase of surface heat transfer, and on the other hand, it reduces the temperature gradients and the intensity of heat transfer. The last control factor to be assessed is the number of undulations N of the wavy tube. The thermal and flow patterns of the PCM in the enclosure are shown in Figures 11 and 12 for two values of N, N = 4 (the optimal value) and N = 1, and for the highest wave amplitude of 3 mm. The flow streamlines mimic the wall profile near the hot tube, as the number of oscillations of the streamlines is equal to that of the tube. The volume of melted PCM, as well as the temperature in different zones of the cavity, are larger for N = 4 compared to N = 1, which can be attributed to the increase of surface contact between the PCM and the hot tube when the same amplitude is considered. This results in an increase of the MVF and the ES when N is raised, as observed in Figure 13. However, the presence of undulation entraps the fluid flow at the tube side and reduces the heat transfer. Thus, on the one hand, the presence of N contributes to the increase of surface heat transfer, and on the other hand, it reduces the temperature gradients and the intensity of heat transfer.

Conclusions
In the present numerical study, an LHTES unit's design based on a wavy tube and NePCM was optimized, based on Taguchi's approach. The control factors considered were the amplitude A of the wavy tube and its number of undulations N, and the volume fraction of the nanoparticles ω na .

•
Each control factor was assigned four different values corresponding to four levels. The L16 orthogonal array method was used to reduce the number of required tests. The time required to melt the PCM fully was selected as the optimization result, so the design was optimized based on: the lower the full melting time, the better the result. • It was found that the optimal values of the control factors are: ω na = 0.06 and a plane tube with no undulation (A = 0 and N = not applicable). Conducting further simulations with values around the optimal ones confirmed that the melting time is the lowest when the optimal values are selected. Plotting the isotherms and the streamlines for various values of the control factors showed that reducing the amplitude improves heat transfer in the bottom region and enhances PCM melting and energy storage.

•
Increasing the nanoparticle fraction elevates the thermal conductivity of the PCM and enhances heat transfer, resulting in an accelerated melting and high stored energy. Finally, increasing the number of undulations of the wavy tube increases the heated contact surface with the PCM, which leads to a faster PCM melting and improves stored energy.