The Effect of Variable-Length Fins and Different High Thermal Conductivity Nanoparticles in the Performance of the Energy Storage Unit Containing Bio-Based Phase Change Substance

Abstract

Thermal Energy Storage (TES) is a key feature in the sizing of thermal systems and energy management. The Phase Change Material (PCM) can store a huge amount of heat in the form of latent heat. However, a good design of the TES unit is required to absorb thermal energy and charge quickly. In the present study, a combination of optimum fin design and nanoadditives are used to design a shell and tube shape TES unit. The Taguchi optimization method is employed to maximize the melting rate by optimizing the arrangement shape of fins and the type and the volume fractions of nanoparticles. The results showed that long fins should be mounted at the bottom and short fins at the top, so that the PCM melts down at the bottom quickly, and consequently, a natural convection circulation occurs at the bottom and advances in the solid PCM. The short fins at the top allow a good natural convection circulation at the top. An increase in the volume fraction of nanoparticles increases the melting rate. An optimum design shows a 20% more melting rate compared to a poor design.

1. Introduction

The design of Thermal Energy Storage (TES) units has gained considerable attention during the past few years since they can hold a significant amount of heat in a compact space [1]. The TES units significantly contribute to the thermal management of devices and heat loading shifts. For example, a TES unit can be integrated into a cooling/heating system subject to fluctuations in thermal loads, so that the TES could absorb the fluctuations, and the thermal system can be designed for the average thermal load. Moreover, the TES systems can shift the thermal loads between hours, days, or even seasons [2]. The TES heatsinks are also a novel type of heatsinks that can effectively damp transient heat loads in electronic components [3,4]. PCMs have also been integrated with bricks [5] to improve the building material’s heat transfer characteristics.

The Phase Change Materials (PCMs) are a typical form of thermal energy storage medium, which could undergo a melting/solidification process and store/release thermal heat in the form of latent heat. PCMs are an excellent candidate for compact thermal energy storage; however, their low thermal conductivity results in poor heat transfer properties [6]. A poor heat transfer rate leads to a slow response time of TES systems, and a long time is required for charging or discharging of thermal energy. Therefore, thermal enhancement approaches are required to improve the heat transfer rate of PCMs, and consequently, the response time of TES. Several approaches such as using metal foams [7], thermal conductivity enhancement of PCMs by nanoparticles [6], and extended surfaces and fins [8] have been proposed to improve the heat transfer of TES units.

The fins have shown effective improvement in the charging/discharging time of TES units since they distribute the heat inside the PCM. Fins increase the contact area of PCM and hot/cold sources and lead to heat transfer improvement. However, the mass of fins does not contribute to latent heat energy storage, and thus, using metallic fins increases the size and weight of a TES for a fixed capacity of thermal energy storage. Moreover, typically the fins are made of copper or aluminum, which are expensive metals. Thus, the mass of utilized fins should be minimized. Another important point in using fins is the fact that they act as obstacles and could suppress the natural convection flows in a molten region of a TES capacity. Considering these points, many researchers tried to investigate the impact of using fins on the heat transfer behavior of TES units.

Borhani et al. [9] analyzed the melting heat transfer of RT50 in a shell-and-tube shape TES unit. The space between the shell and the tube was filled with a PCM. A hot fluid was flowing in the tube while the surface of the tube was supported with a spiral-shaped fin extended inside the PCM. The authors studied the impact of a geometrical design of fins, such as the fin’s pitch and thickness, on the charging process. The outcomes showed that increasing the fin pitch from 10 mm to 20 mm, decreases the melting time by 35%. A good design of thickness and pitch enhanced the melting performance by 56%. Deng et al. [10] examined the melting heat transfer of lauric acid PCM in a shell and tube heat exchanger TES unit. They mounted two fins at the tube side into the PCM and investigated the impact of placement location (angle) or fins. They found that there is an optimum angle for the placement of fins that could increase the melting performance by 66.7%.

The combination of metal foams and PCMs [11] is another approach that has been investigated recently. For example, Talebizadehsardari et al. [12] examined the discharging of an air heater filled with PCM-metal foam composite. Authors found that a well-designed geometry of heat transfer passages can change the discharging time from 13.6 h to 6 h. The impact of the storage size [13] on the thermal energy storage and discharging behavior [14] in domestic applications have also been examined.

Some researchers tried to promote the heat transfer of PCMs by adding nanoadditives to phase change materials and produce nano-enhanced phase change materials (NePCMs) with improved thermal conduct. Sheikholeslami et al. [15] employed a Snowflake-shape fin configuration to improve the solidification of a PCM in a TES unit and reported that the presence of the Snowflake fins could notably increase the solidification rate without imposing a significant influence on the thermal energy storage capacity of the unit. Ma et al. [16] employed a NePCM to improve the solidification rate of PCM in a wavy thermal energy storage duct. The duct walls were filled with NePCM while the cold air was blowing into the wavy surface of the channel. The authors reported that the increase of the amplitude of channel waves reduces the discharging time of the TES unit. Moreover, the presence of nanoadditives were pronounced when the channel surface was wavy. The impact of using nanoadditives for heat transfer in channel shape TES units [17] has also been examined. Very recently, Irwan et al. [18] reviewed the phase change process of NePCMs.

The literature review shows that either using fins or nanoadditives could improve the charging/discharging time of a TES unit. The combination of these two approaches could also be effective. It is clear that the presence of nanoparticles changes the thermophysical properties of PCM while the fins could suppress the convective flows in molten areas of TES. However, there is not any research that systematically evaluates a design of TES units to fully benefit from the advantage of nanoadditives and the geomatical shape of fins. Thus, the present research aims to address the impact of using fins and various types of NePCMs to minimize the charging time of a TES unit. The Taguchi optimization approach would be employed to find the optimum type of nanoparticles and fins arrangements.

2. Mathematical Model

2.1. Model Description

Figure 1a–c depicts 3D geometry shell and tube as a latent heat energy storage unit with annular fins (rectangular 2D cross-section) in different configurations. In Figure 1a, the length of the fins is taken to be the same. However, in Figure 1b,c, the length of the located fins along the heat transfer flow direction can be descending and increasing. The ratio of the length of the lowest fin to the base fin is AR. The base fin refers to the fin located in the configuration of Figure 1a. It is worth noting that the used fins’ volume in all structures is constant and considered as the designing constraint. The geometric characteristics of the unit are as follows: The internal tube, with a radius of 5 mm and ignorable thickness, is made of copper. The thin fins made of copper have a thickness of 1 mm and variable lengths. Geometrical characteristics of the unit can be seen in Figure 1d. The heat transfer fluid (HTF) flowing inside the internal tube is water, and the base PCM of the composite is coconut oil. The dispersed nanoparticles as the high thermal conductivity materials can be Al₂O₃, graphene oxide (GO), Cu, or Ag. The specifications of the coconut oil, the copper fin, the HTF, and the nanoparticles can be found in

Table 1. Thermophysical properties of the coconut oil and the nanoadditives [19,20].

Properties	Copper Fin	Water	Coconut Oil (Measured)
			Solid (15 °C)	Liquid (32 °C)
$\rho$ (kg m−3)	8900	993.73	920	914 ± 0.11%
Cp (J kg K−1)	386	4178	3750	2010 ± 0.2%
η (N s m−2)	-	0.705 × 10−3	-	0.0326 ± 3%
k (W m−1 K−1)	401	0.623	0.228	0.166 ± 1.2%
hf (kJ kg−1)	-			103 ± 1%
Pr	-			394.73 ± 3.2%
Tme (ΔTme)	-		24 °C (±1 °C)

and

Table 2. Nanoparticles dispersed in the phase change material.

Properties of Nanoadditives	Al2O3	GO	Ag	Cu
${\rho }_{PCM}$ (kg m−3)	3600	1800	10,500	8960
Cp (J kg−1 K−1)	765	717	235	385
k (W m−1 K−1)	36	5000	429	400
β (K−1)	7.8 × 10−6	28.4 × 10−5	18.9 × 10−6	16.7 × 10−6

.

2.2. 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:

Mass conservation

$$\frac{q_r}{r}+\frac{\partial q_r}{\partial r}+\frac{\partial q_z}{\partial z}=0$$

(1)

Momentum equations

$$\begin{array}{l}{\rho }_{NP,l}\left(\frac{\partial q_r}{\partial t}+q_r\frac{\partial q_r}{\partial r}+q_z\frac{\partial q_r}{\partial z}\right)=-\frac{\partial p}{\partial r}+\frac{{\eta }_{NP,l}}{r}\frac{\partial q_r}{\partial r}+{\eta }_{NP,l}\left(\frac{{\partial }^2{q}_r}{\partial r^2}+\frac{{\partial }^2{q}_r}{\partial z^2}\right),\\ -\frac{{\eta }_{NP,l}{q}_r}{r^2}+\frac{{\omega }^2\left(T\right)-2\omega \left(T\right)+1}{\lambda +{\omega }^3\left(T\right)}{A}_{mush}{q}_r\end{array}$$

(2a)

$$\begin{array}{l}{\rho }_{NP,l}\left(\frac{\partial q_z}{\partial t}+q_r\frac{\partial q_z}{\partial r}+q_z\frac{\partial q_z}{\partial z}\right)=-\frac{\partial p}{\partial z}+\frac{{\eta }_{NP,l}}{r}\frac{\partial q_z}{\partial z}+{\eta }_{NP,l}\left(\frac{{\partial }^2{q}_z}{\partial r^2}+\frac{{\partial }^2{q}_z}{\partial z^2}\right),\\ +g{\rho }_{NP,l}{\beta }_{NP,l}\left(T-T_{me}\right)+\frac{{\omega }^2\left(T\right)-2\omega \left(T\right)+1}{\lambda +{\omega }^3\left(T\right)}{A}_{mush}{q}_z\end{array}$$

(2b)

The 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⁵ kg/m³s)) and small (10⁻³) values, respectively. The source term $\left({\omega }^2\left(T\right)-2\omega \left(T\right)+1/\lambda +{\omega }^3\left(T\right)\right)A_{mush}q$ simulates a resistance force similar to that of a Darcy porous space with variable permeability. The resistance force is a function of the local liquid fraction (ω(T)) and fluid velocity. This term approaches zero in the liquid region and disappears. In a solid region where the liquid fraction is zero, ω(T) = 0, this term approaches (λ × A_mush) × q, which leads to considerable resistance force agents in the material motion (q) and forces the velocities to zero in solid PCM. This source term changes the force smoothly from zero in liquid to a considerable value in solid. The continuous and smoothness of the force variation is important for numerical convergence of solver. Function ω(T) of the above equations is described as the following:

$$\omega \left(T\right)=\{\begin{array}{ll}0& T<T_{me}-\Delta T_{me}/2\\ \frac{T-T_{me}}{\Delta T_{me}}+\frac{1}{2}& T_{me}-\Delta T_{me}/2<T<T_{me}+\Delta T_{me}/2\\ 1& T>T_{me}+\Delta T_{me}/2\end{array},$$

(3)

Energy conservation

$$\begin{array}{l}{\left(\rho C_p\right)}_{NP}\frac{\partial T}{\partial t}+(\rho C_p){}_{NP,l}\left(q_r\frac{\partial T}{\partial r}+q_z\frac{\partial T}{\partial z}\right)=\frac{1}{r}\frac{\partial }{\partial r}\left(k_{NP}r\frac{\partial T}{\partial r}\right)\\ +\frac{\partial }{\partial z}\left(k_{NP}\frac{\partial T}{\partial z}\right)-\left(1-V{F}_{na}\right){\rho }_{PCM,l}{h}_{f,PCM}\frac{\partial \omega \left(T\right)}{\partial t}\end{array}$$

(4)

in which,

$${\left(\rho C_p\right)}_{NP}=\omega \left(T\right)\left[{\left(\rho C_p\right)}_{NP,l}-{\left(\rho C_p\right)}_{NP,s}\right]+{\left(\rho C_p\right)}_{NP,s},$$

(5a)

$$k_{NP}=\omega \left(T\right)\left(k_{NP,l}-k_{NP,s}\right)+k_{NP,s},$$

(5b)

Energy conservation for solid fins

$${\left(\rho C_p\right)}_s\frac{\partial T}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(k_{s}r\frac{\partial T}{\partial r}\right)+\frac{\partial }{\partial z}\left(k_s\frac{\partial T}{\partial z}\right),$$

(6)

2.3. The Thermophysical Properties of Nano PCM

The following relation was used to evaluate the NePCM density:

$${\rho }_{NP}={\rho }_{PCM}+{\nu }_{na}\left({\rho }_{na}-{\rho }_{PCM}\right),$$

(7)

where in the above equation, the concentration of nanoparticles is denoted by ${\nu }_{na}$. The effective dynamic viscosity, thermal conductivity, and thermal expansion of the NePCM are evaluated as:

$${\eta }_{NP,l}={\eta }_{PCM,l}{\left(1-{\nu }_{na}\right)}^{-2.5},$$

(8)

$$\frac{k_{NP,i}}{k_{PCM,i}}=\frac{\left(k_{na}+2k_{PCM,i}\right)-2{\nu }_{na}\left(k_{PCM,i}-k_{na}\right)}{\left(k_{na}+2k_{PCM,i}\right)+{\nu }_{na}\left(k_{PCM,i}-k_{na}\right)},$$

(9)

$${\rho }_{NP,l}{\beta }_{NP,l}={\rho }_{PCM,l}{\beta }_{PCM,l}+{\nu }_{na}\left({\rho }_{na}{\beta }_{na}-{\rho }_{PCM,l}{\beta }_{PCM,l}\right),$$

(10)

where in the above index, i shows the phase of NePCM, which could be solid or liquid. Besides, the effective heat capacity is a function of volume fractions of nanoparticles and was estimated as:

$${\left(\rho C_p\right)}_{NP}={\rho }_{PCM}{C}_{p,}{}_{PCM}+{\nu }_{na}\left({\rho }_{na}{C}_p{}_{,na}-{\rho }_{PCM}{C}_{p,}{}_{PCM}\right),$$

(11a)

$${\rho }_{PCM}{C}_{p,}{}_{PCM}\left(T\right)=\omega \left(T\right)\left[{\rho }_{PCM,l}{C}_p{}_{,PCM,l}-{\rho }_{PCM,s}{C}_p{}_{,PCM,s}\right]+{\rho }_{PCM,s}{C}_p{}_{,PCM,s},$$

(11b)

2.4. 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:

Conservation mass

$$\frac{q_r}{r}+\frac{\partial q_r}{\partial r}+\frac{\partial q_z}{\partial z}=0,$$

(12)

Momentum equations

$${\rho }_{water}\left(\frac{\partial q_r}{\partial t}+q_r\frac{\partial q_r}{\partial r}+q_z\frac{\partial q_r}{\partial z}\right)=-\frac{\partial p}{\partial r}+\frac{{\eta }_{water}}{r}\frac{\partial q_r}{\partial r}+{\eta }_{water}\left(\frac{{\partial }^2{q}_r}{\partial r^2}+\frac{{\partial }^2{q}_r}{\partial z^2}\right)-\frac{{\eta }_{water}{q}_r}{r^2},$$

(13a)

$${\rho }_{water}\left(\frac{\partial q_z}{\partial t}+q_r\frac{\partial q_z}{\partial r}+q_z\frac{\partial q_z}{\partial z}\right)=-\frac{\partial p}{\partial z}+\frac{{\eta }_{water}}{r}\frac{\partial q_z}{\partial r}\frac{\partial }{\partial r}+{\eta }_{water}\left(\frac{{\partial }^2{q}_z}{\partial r^2}+\frac{{\partial }^2{q}_z}{\partial z^2}\right),$$

(13b)

Energy conservation

$${(\rho C_p)}_{water}\left(\frac{\partial T}{\partial t}+q_r\frac{\partial T}{\partial r}+q_z\frac{\partial T}{\partial z}\right)=\frac{1}{r}\frac{\partial }{\partial r}\left(k_{water}r\frac{\partial T}{\partial r}\right)+\frac{\partial }{\partial z}\left(k_{water}\frac{\partial T}{\partial z}\right),$$

(14)

2.5. Controlling Boundary and Initial Conditions

The boundary conditions at the interface of the wavy wall and NePCM are the continuities of the temperature and heat flux as the following:

$${T|}_{HTF}={T|}_{NP},{k_{HTF}\frac{\partial T}{\partial n}|}_{HTF}{}_{}={k_{NP}\frac{\partial T}{\partial n}|}_{NP},$$

(15a)

The HTF is injected with the temperature of T_in and v_in.

$${T|}_{water}=T_{in},{q_r|}_{water}=0,{q_z|}_{water}=q_{z,in},$$

(15b)

At the outlet of the HTF tube, the developed convective flow is considered to be established.

$${q_r|}_{water}=0,{\frac{\partial T}{\partial z}|}_{water}{}_{}={\frac{\partial q_z}{\partial z}|}_{water}=0,$$

(15c)

Due to the symmetrical distribution of the tube bundle, the boundary conditions of the computational domain are:

$${q_r|}_{NP}={q_z|}_{NP}=0,{\frac{\partial T}{\partial r}|}_{NP}{}_{}=0,$$

(15d)

At the upper and lower surfaces of the tube bundle, the boundary conditions can be defined as the following:

$${q_r|}_{NP}={q_z|}_{NP}=0,{\frac{\partial T}{\partial z}|}_{NP}{}_{}=0,$$

(15e)

The initial condition for the NePCM domain can be expressed as the following:

$${q_r|}_{NP}={q_z|}_{NP}=0,{T|}_{NP}{}_{}=T_{initial},$$

(15f)

2.6. Total Energy Stored in the Unit (ES) and Melting Rate (MVF)

The cumulative heat transferred from the HTF to the NePCM is exactly the stored energy in the unit and can be calculated by integrating the instantaneous heat transfer rate from the wavy wall:

$$TotalEnergy={\displaystyle \underset{A}{\int }{\left(\rho C_p\right)}_{NP}\left(T-T_{in}\right)}dA+{\displaystyle \underset{A}{\int }\left(1-V{F}_{na}\right){\rho }_{PCM,l}{h}_{f,PCM}dA},$$

(16)

the melting volume fraction is evaluated as:

$$MVF\left(t\right)=\frac{{\displaystyle {\int }_A^{}\omega \left(T\right)dA}}{{\displaystyle {\int }_A^{}dA}},$$

(17)

To further discuss the temperature distribution, the temperature uniformity index σ°(t) was proposed, as described below:

$$\sigma \left(t\right)={\left[{\displaystyle \underset{\forall }{\int }{\left(T_{av}-T\left(r,\theta ,z,t\right)\right)}^{2}d\forall }\right]}^{0.5}{\left({\displaystyle \underset{\forall }{\int }d\forall }\right)}^{-1},$$

(18)

$$T_{av}=\left({\displaystyle \underset{\forall }{\int }Td\forall }\right){\left({\displaystyle \underset{\forall }{\int }d\forall }\right)}^{-1},$$

(19)

3. Numerical Approach, Grid Dependency and Validations

In this investigation, controlling equations for the phase change flow are solved by employing the finite element numerical approach implemented in the Galerkin technique, coded as user-defined functions. The transitional forms of the governing equations, called weak forms, are established in the phase change material zone, the fins zone, and the HTF zone. To approximate the field variables, Lagrange elements of various orders are employed. To iteratively solve the residual equations, controlling the convergence, the PARallel DIrect SOlver (PARDISO) [21,22] is employed. The applied Newtonian damping factor for the solver is 0.8 and the relative tolerance is 10⁻⁴. Backward differentiation formula (BDF) is utilized to handle the time dependent computations. More details about the utilized method and solution approach could be found in [23,24].

The computational domains defined for HTF and fins are meshed by structural rectangular grids. However, un-structural triangular grids are employed to mesh the NePCM domain. A useful tactic to reach an accurate solution with a low computational cost is the use of an optimum grid having a high-quality. Various meshes with different grid sizes are applied to check the impact of grid sizes on the simulated results. Several meshes with different grid sizes are applied to check the impact of grid sizes on the simulated results.

Table 3. Different mesh sizes.

Case Study	C1	C2	C3	C4	C5
No. of Elements	15,200	22,440	35,550	51,660	81,600
Computation time (min)	3 h 57	5 h 26	8 h 57	18 h 55	19 h 17
Temperature point#1(K) at 1500 s	295.9986	295.9984	295.9991	295.9985	295.999
Temperature point#2 (K) at 1500 s	300.1385	300.1381	300.1351	300.1356	300.1344
MVF at 1500 s	0.83112	0.83102	0.83101	0.83097	0.83096

lists the grids with the different number of elements used for meshing. Figure 2a,b illustrates the dependency of the MVF and temperature calculated at P₅ to the mesh size. As these figures depict, no considerable changes in the MVF and temperature are observed when changing the mesh sizes. Although using a finer mesh could increase the computational cost, the solver was sensitive to the mesh size, and for a coarse mesh, one should take small steps to keep the convergence of numerical solutions. Indeed, the melting takes place in a narrow space and the distance between mesh elements should be adequately small so the elements should capture the changes. Hence, the mesh of C4 was computed smoothly, and therefore, the grid with 51,660 elements is chosen as the appropriate grid for the simulations.

To assess the reliability of the numerical algorithm, the numerical and experimental studies discussed in [25] are employed. The current results are first compared with those reported in experimental-numerical work performed by Kumar et al. [25]. In their study, the analysis of the melting flow of lead inside a square enclosure, experimentally and numerically, was conducted. Figure 3 illustrates the progressive melting fronts over time for the conducted work by Kumar et al. [25] and the present model.

In the second comparative analysis, the temperature distribution resulted from the current simulation is evaluated by employing that reported in [26] and plotted in Figure 4. This figure shows the isolines of natural convection in a portioned cavity heated from the left and cooled from right. The top and bottom are well insulated. The regular enclosure containing two vertical fins was filled by air.

Figure 5 shows the melting fraction for in a cavity filled with 1% wt. graphene nanoplatelets dispersed in 1-tetradecanol as experimented in [27]. The cavity was rectangular with a size of (25 mm) × width (20 mm), which led to an aspect ratio of 0.8. The NePCM was at an initial super cold temperature of 36 °C. During the melting, the heated wall was kept at an isothermal temperature of 47 °C while the fusion temperature was 37 °C. The melting fraction as a function of time is compared with the numerical results of the present study. The three comparative analyses discussed above confirm the reliability of the employed model.

4. Results and Discussion

The design parameters of the TES unit are the volume concentration of nanoparticles (0 < υ_na < 4.5%), the aspect ratio of the PCM enclosure (0.1 < AR < 1.45), and the type of nanoparticles. The Taguchi optimization method is adopted to explore the combination of possible designs and find the best values of these design parameters to maximize the melting rate. Following the Taguchi method, each design parameter should be divided into levels. Here, four levels are selected for each parameter. The selected levels are depicted in

Table 4. The range and levels of control parameters.

Factors	Description	Level 1	Level 2	Level 3	Level 4
A	υna, volume fraction (%)	0	1.5	3	4.5
B	AR, radius ratio	0.1	0.55	1	1.45
C	type of nanoparticles	Al2O3	GO	Ag	Cu

. Four types of nanoparticles are also selected, which are Al₂O₃, GO, Ag, and Cu.

The Taguchi method uses orthogonal tables to efficiently explore the design space. Since there are three design parameters and four levels, the possible design combinations are 4³. Here the standard orthogonal table of L16 was selected, which is summarized in

Table 5. Taguchi table L16 for orthogonal search of design space.

Case	A	B	C	Value at 19,800 s			S/N Ratio
	υna	AR	Nano	MVF	Total Energy (J)	σ
1	0.0	0.10	Al2O3	0.782	165,239	0.298	−2.13433
2	0.0	0.55	GO	0.854	177,186	0.271	−1.37403
3	0.0	1.00	Ag	0.917	186,140	0.234	−0.74990
4	0.0	1.45	Cu	0.969	189,298	0.198	−0.27462
5	1.5	0.10	GO	0.796	166,329	0.298	−1.98437
6	1.5	0.55	Al2O3	0.866	178,470	0.269	−1.25042
7	1.5	1.00	Cu	0.927	187,426	0.230	−0.65516
8	1.5	1.45	Ag	0.974	190,179	0.186	−0.22818
9	3.0	0.10	Ag	0.808	167,881	0.297	−1.84641
10	3.0	0.55	4 Cu	0.878	180,098	0.266	−1.13340
11	3.0	1.00	Al2O3	0.937	188,172	0.225	−0.56629
12	3.0	1.45	GO	0.979	190,429	0.176	−0.18192
13	4.5	0.10	Cu	0.820	169,795	0.294	−1.71961
14	4.5	0.55	Ag	0.890	180,740	0.263	−1.00946
15	4.5	1.00	GO	0.948	188,335	0.219	−0.46525
16	4.5	1.45	Al2O3	0.984	191,622	0.169	−0.13927

. The volume concentration for cases 1–4 of Table 5 are zero, and thus, there is no difference between the type of nanoparticles. However, the aspect ratio is changing and thus these cases are not identical.

The simulations were executed for each case of Table 5, and the melting rate, stored energy, and temperature uniformity of the PCM enclosure were computed. The melting time of 5.5 h (19,800 s) was adopted for the design target. This is a typical time for charging a TES unit.

The computed values of the molten portion (MVF), total stored energy (ES), and temperature uniformity were written in Table 5. The molten portion (MVF) was adopted as the target parameter with the aim of “the larger, the better”. Using the Taguchi method, the values of signal to noise (S/N) were computed and inserted in Table 5. A large value of S/N shows the advantage of a case for the selected target (higher MVF). Now, the computed S/N ratios were utilized to plot S/N ratios for each design parameter, as illustrated in Figure 6.

Following the Taguchi method, a level corresponding to the highest value of S/N for each design parameter should be selected as the optimum level. Thus, υ_na = 0.045, AR = 1.45, and Cu nanoparticles were selected as the optimum design case.

Table 6. The optimum values of the controlling parameters.

Optimum Factors			Optimum Melting at 19,800 s
υna	AR	Nano	Taguchi Prediction	Tested Case
0.045	1.45	Cu	0.995	0.984

shows the estimated MVF at 19,800 s, which is 0.995. The simulations were done for the same case, and the actual MVF was computed as 0.984, which is close to the estimated value. This case is very close to case 16 of Table 5 for Al₂O₃. However, the Cu nanoparticles could be more compatible with the copper material of fins, and thus, the case of Cu nanoparticles was selected as the optimum nanoparticles. The optimum case shows a 20% more melting rate than a poorly designed case of 1 (from Table 5).

The S/N ratio denotes a design advantage for a test case. The larger S/N ratio, the more influential the parameter’s level. Moreover, the larger the variation of S/N, the more effective the parameter. Figure 6 shows that the variation of S/N is highest for aspect ratio (AR) and lowest for the type of nanoparticles. The volume fraction of nanoparticles is the second most important, influential design parameter. It is clear that the increase of volume fraction of nanoparticles and aspect ratio increases the MVF.

A further investigation of the impact of variation of volume fractions of Cu nanoparticles and aspect ratio on the melting behavior was done. The investigation is designed to explore the melting behavior around the optimum design case. Thus, six cases of

Table 7. Further investigation on the impact of design parameters around the optimum point for Cu nanoparticles.

Case	Parameter Investigation	A	C	Value at 19,800 s
		υna	AR	MVF	Total Energy (J)	σ1/2
1	υna	0	1.45	0.96889	189,304	0.1976
2		0.015	1.45	0.97399	190,424	0.18642
3		0.03	1.45	0.97876	191,312	0.17691
4	AR	0.045	0.1	0.82042	169,804	0.29451
5		0.045	0.55	0.88912	181,427	0.26202
6		0.045	1	0.94594	189,817	0.21876
7	Optimum	0.045	1.45	0.98399	192,316	0.16868

were adopted, and the melting process for these cases was simulated.

Figure 7 shows the temperature distribution in the enclosure for various volume fractions of nanoparticles and different melting times when AR = 1.45. At a short time of 1200 s, the tube’s solid wall is at a hot temperature, and almost the entire enclosure is at a cold temperature. The heat transfer between HTF fluid and PCM in the enclosure develops mostly at the bottom parts of PCM. In the beginning of heat transfer, the PCM is in a super-cold solid state. Thus, there is no natural convection circulation in the PCM domain. The HTF fluid enters the TES unit as the bottom, and hence, the most robust heat transfer rate occurs at the bottom. The PCM around the fins is hot, which shows the impact of fins at this stage. Moreover, since AR = 1.45, the fins at the bottom are much longer than the top. At the time step of 1200 s, the impact of variation of υ_na on the temperature field is not obvious. As time passes, the melting advances, and the temperature gradients extend to an un-melted area. The impact of the presence of nanoparticles becomes obvious after long periods of time, where the melting advances.

Figure 8 shows the streamlines and melting interface in the PCM domain. The melting starts around the fins and advances inside the enclosure during the time. A comparison between the parents on melting interface and solid regions for various concentrations of nanoparticles shows that the increase of concentration increases the melting rate and reduces the remaining solid regions. Figure 9a,b depicts the MVF and total stored energy during the melting process. Figure 9a shows that the higher the concentration of nanoparticles, the larger the melting rate. An increase in nanoparticle concentration also decreases the charging time. These findings are in good agreement with the melting shapes observed in the contours of Figure 8. Figure 9b also indicates that the rise of nanoparticle concentration increases the rate of stored energy. The presence of nanoparticles mostly influences the stored energy in the middle of the charging process.

Figure 10 displays the temperature evolution of monitoring points in the enclosure during the melting process. This figure shows that there are three main stages for temperature evolution of monitor points. First, there is a semi-linear increase of temperature, which is related to the pure conduction regime in the supper cold PCM. Then, there is a semi-flattened temperature change in which the temperature stays constant at a temperature about the fusion temperature. This is where the phase change occurs, and the energy melts the PCM instead of raising its temperature. After this stage, there is a semi-linear increase of temperature, which corresponds to the temperature variation in the molten PCM and in the presence of natural convection effects. Figure 10 shows the impact of variation of Cu nanoparticles on the temperature profiles. As seen, at the first stage, which is pure conduction, the increase of nanoparticle concentration slightly increases the temperature of monitor points. The presence of the nanoparticles increases the thermal conductivity of the NePCM, and thus the heat could be more easily transfer from the HTF tube inside the PCM. Thus, at a constant time, the temperature of the monitor points is higher for a larger volume fraction of nanoparticles. The presence of nanoparticles does not change the magnitude of the second part (the melting part) of the profiles since the fusion temperature does not change. However, the increase of concentration of nanoparticles accelerates the beginning of the third stage (heat transfer in molten PCM). Two reasons could contribute positively to this observation. First, an increase in nanoparticles concentration improves the heat transfer rate, and thus the melting finishes faster. The second reason is the fact that the presence of nanoparticles slightly reduces the latent heat of phase change, which consequently leads to a faster phase change.

Figure 11 shows the temperature uniformity index (σ) in the PCM region during the phase change for various concentrations of nanoparticles. The increase of concentration of nanoparticles induces minimal impact on the temperature uniformity of TES in the early stages of melting. Indeed, in Figure 10, it was observed that the variation of nanoparticle concentration slightly changes the temperature profiles in the solid PCM and during the melting stage. Thus, this observation is in agreement with Figure 9. There is a notable difference in the temperature uniformity of TES at the final stage of melting. This is because at the end of the melting process, where there is no or minimal solid PCM remaining, the overall temperature of the enclosure tends to the temperature of HTF fluid. The full charging occurs faster for a high concentration of nanoparticles, and thus, σ also drops quickly.

Figure 12 shows the isotherm contours of the TES unit for various aspect ratios (AR). When the AR is smaller than unity, the bottom fins are short, and the top fins are long. In the case of AR = 1, the length of the fins is equal. For the case of AR > 1, the bottom fins are longer than the top fins. The total mass of the fins is constant in all cases, and the only difference is the variation of the geometrical shape and arrangement. Figure 12 shows the significant impact of AR on the temperature distribution in the enclosure. Figure 13 illustrates the melting contours and streamlines in the enclosure for various aspect ratios. When AR = 0.1, it is clear that the melting advances to the right from the top, where the fins have the highest length. Thus, the top region is fully melted at the final stage of melting while the bottom region is at a solid-state. The hot molten PCM moves toward the top, where there is no PCM to be melted. However, when AR = 1.45, the long bottom fins melt a significant amount of PCM at the bottom. Then, the natural convection heat transfer occurs. The hot molten PCM moves toward the top region where there is plenty of cold PCM and contributes to the melting of the top regions of the enclosure where the fins were short. Thus, the natural convection notably aids the melting of PCM in the case of AR = 1.45.

Figure 14 illustrates the melting profiles and stored energy for various values of AR. As seen, the increase in AR increases both MVF and stored energy. The depicted melting behavior of Figure 13 confirms the behavior of this figure very well. Since the concentration of nanoparticles is fixed as υ_na = 0.045, the increase of MVF could directly increase the total stored energy. These could be slight fluctuations that are due to the sensible part of stored energy, and their impact on toral stored energy is minimal.

Figure 15 shows the temperature profiles of monitor points for various values of AR. The bottom point is under significant influence of AR, since the bottom is a conduction-dominant region, and the variation of the length of the fins directly affects the conduction heat transfer in this part of the enclosure. The temperature history of P3 shows that the melting behavior of this point is almost independent of AR. This is because the variation of changes in the fin’s lengths next to this point is minimal. The temperature of P5 in the liquid region is the lowest for case AR = 1.45 because the melting rate for this case is maximum, and the transferred energy from the HTF fluid mostly consumes by the latent heat. P8 shows a smooth transition between stages since this point is at the top, and regardless of the length of the fins, there is a slight natural convection effect at this point, and the PCM will melt quickly.

Figure 16 displays the temperature non-uniformity (σ) for various values of AR. The variation of fins shape (AR) induces minimal impact of σ at early times. This is because all PCM is solid at first, and the heat transfer is conduction dominant. Thus, the location of fins is not important in the general behavior of heat transfer. However, in the middle and final stages of charging, the natural convection flows occur, and the arrangement of fins lead to significant differences. At the final stage of charging, the TES unit with AR = 1.45 charges faster, and thus, the temperature of molten PCM quickly approaches the HTF temperature, and σ drops quickly.

5. Conclusions

A shell-and-tube shape thermal energy storage unit was optimized for the shape of fins, type of nanoparticles, and volume fraction of nanoparticles. The amount of the fins’ mass was considered fixed, but the length of the fins could change from bottom to top. The natural convection effects in the molten region were taken into account. The finite element method was applied to integrate the continuity, momentum, and phase change convective heat transfer partial differential equations in an axis-symmetric cylindrical coordinate system. The mesh study and verifications were performed, and the model was found in fair agreement with experimental and numerical literature works. The Taguchi optimization method was applied to maximize the melting rate in the enclosure. The melting process was investigated in detail through melting and temperature distribution contours, streamlines, and graphs. The following key points were found:

• The uniform shape of the fins was important since the fins could have several important impacts on the melting process. The presence of the fins supports the pure conduction regions and improves the conduction heat transfer. However, the fins also resist the natural convection circulation and tend to suppress the convection flows. Moreover, the natural convection flows tend to melt top regions rather than the bottom. Thus, the fins should help with melting at the bottom. Thus, the long fins at the bottom and short fins at the top were advantageous. The long fins at the bottom melt the PCM by conduction mechanism effectively. The short fins at the top also allow circulation of molten PCM and the top regions.
• The types of nanoparticles do not induce a significant impact on the melting rate and stored energy. The Cu and Al₂O₃ particles show almost the same optimum performance. However, the fins and enclosure walls are made of copper, and hence, here, the copper nanoparticles were adopted as the optimum material.
• The increase of the volume fraction of nanoparticles could improve the heat transfer and accelerate the charging process.
• The present investigation results showed that the fins’ geometrical arrangement was the most influential design parameters to increase the melting rate. After that, the volume fraction of nanoparticles was the second important design parameter. The type of nanoparticles was the least effective parameter. The optimum design could increase the melting rate by 20%.

It should be noted that the synthesis and stabilization of NePCMs add some degree of complexity to a TES design. Hence, the minor advantage of using NePCMs could be investigated regarding the TES unit’s cost and complexity. Moreover, the impact of the shape, average size of the nanoparticles, and additives on the unit’s thermal performance were neglected. Such effects could better highlight the difference between different types of nanoparticles and could be subject to future studies.

The geometrical arrangement of the fins was the most influential design parameters to increase the melting rate. After that, the volume fraction of nanoparticles was the second important design parameter. The type of nanoparticles was the least effective parameter. The optimum design could increase the melting rate by 20%.