Optimization of Nano-Additive Characteristics to Improve the Efﬁciency of a Shell and Tube Thermal Energy Storage System Using a Hybrid Procedure: DOE, ANN, MCDM, MOO, and CFD Modeling

: Using nano-enhanced phase change material (NePCM) rather than pure PCM signiﬁcantly affects the melting/solidiﬁcation duration and the stored energy, which are two critical design parameters for latent heat thermal energy storage (LHTES) systems. The present article employs a hybrid procedure based on the design of experiments (DOE), computational ﬂuid dynamics (CFD), artiﬁcial neural networks (ANNs), multi-objective optimization (MOO), and multi-criteria decision making (MCDM) to optimize the properties of nano-additives dispersed in a shell and tube LHTES system containing parafﬁn wax as a phase change material (PCM). Four important properties of nano-additives were considered as optimization variables: volume fraction and thermophysical properties, precisely, speciﬁc heat, density, and thermal conductivity. The primary objective was to simultaneously reduce the melting duration and increase the total stored energy. To this end, a ﬁve-step hybrid optimization process is presented in this paper. In the ﬁrst step, the DOE technique is used to design the required simulations for the optimal search of the design space. The second step simulates the melting process through a CFD approach. The third step, which utilizes ANNs, presents polynomial models for objective functions in terms of optimization variables. MOO is used in the fourth step to generate a set of optimal Pareto points. Finally, in the ﬁfth step, selected optimal points with various features are provided using various MCDM methods. The results indicate that nearly 97% of the Pareto points in the considered shell and tube LHTES system had a nano-additive thermal conductivity greater than 180 Wm − 1 K − 1 . Furthermore, the density of nano-additives was observed to be greater than 9950 kgm − 3 for approximately 86% of the optimal solutions. Additionally, approximately 95% of optimal points had a nano-additive speciﬁc heat of greater than 795 Jkg − 1 K − 1 . the of conﬁgurations in a and examined in a


The Importance of TES
In recent decades, the environmental challenges associated with energy production using nonrenewable resources combined with the increasing consumption of energy have resulted in a dramatic increase in the use of renewable energy resources. Solar energy is one of the most attractive energy resources for investment in this regard as its annual production capacity is several times that of the world's energy consumption [1,2]. On the other hand, solar energy is a time-dependent resource that cannot operate continuously throughout the day. As a result, thermal energy storage (TES) systems are required to bridge the gap between energy supply and demand and ensure the performance stability of solar energy systems. There are three types of thermal energy storage: thermochemical, latent, and sensible [3]. Latent heat TES has a greater energy storage capacity than sensible heat TES and is more stable than thermochemical TES [4]. Consequently, LHTES systems containing PCM have gained popularity in the field of TES; PCM is used in a wide variety of engineering applications, including energy conservation in buildings [5][6][7], solar collector systems [8][9][10], solar building [11][12][13][14], and photovoltaic systems [15][16][17].

LHTES Units
Recent years have seen a significant increase in numerical and experimental investigations of LHTES units with various geometries. Mahdi et al. [18] numerically simulated the PCM phase change process in a double-pipe helical coiled tube LHTES. They discovered that their proposed LHTES system melts at a rate approximately 25.7% and 60% slower than horizontal and vertical straight double-pipe LHTES systems, respectively. Moreover, their findings indicated that the coil pitch, inlet temperature, and Reynolds number of the heat transfer fluid (HTF) significantly affect PCM melting behavior. Mahdi et al. [19] conducted additional research in which they examined the melting performance of a rectangular LHTES unit numerically and experimentally. Their findings confirmed that partitioning the LHTES unit into eight equal partitions significantly improves natural convection and accelerates the melting process by 65%. Motahar [20] investigated the melting heat transfer and solid-liquid interface tracking in a rectangular unit containing n-octadecane as a PCM. He proposed a multilayer perceptron feedforward neural network for forecasting the Nusselt number and melt fraction during the melting process using Rayleigh, Fourier, and Stefan numbers. The observed results established that natural convection is responsible for a significant portion of the melting process. Mahdi et al. [21] developed a CFD approach to investigate the thermal efficiency of various tube configurations in a shell and tube LHTES unit. They examined three HTF tubes in a variety of configurations within the PCM shell. Their observations revealed that the arc array of tubes at the shell bottom is highly efficient and reduces the melting time by approximately 76% compared to a standard arrangement.
The dispersion of nano-additives in PCM is a widely accepted technique for increasing the efficiency of TES systems. By dispersing high thermal conductivity particles such as metal, metal oxide, and carbon-based nanomaterials into pure PCM, the thermophysical properties of the material are altered, and a range of positive and negative effects on the charging/discharging time and energy storage capability of the material is observed. To this end, Yadav and Sahoo [53] investigated the phase change process in a cylindrical TES system using alumina (Al 2 O 3 )-based NePCM. They considered the effects of total stored energy, the Nusselt number, the Fourier number, and the Rayleigh number. Their findings indicated that NePCM with a 0.3% mass fraction of nano-additives could substantially reduce charging time compared to primary PCM. Chen et al. [54] established an experimental setup to determine the effect of various nano-additives on the performance of solar thermal storage. They dispersed three types of metal and metal oxide nano-additives (Cu, Al 2 O 3 , and TiO 2 ) in PCM at varying volume fractions. Their findings demonstrated that nano-additives could enhance the thermal properties of PCM. Additionally, they stated that PCM containing 0.1% TiO 2 had the lowest viscosity and the best latent heat property. Kumar et al. [55] investigated the effect of CuO nanoparticles on the thermal efficiency of a heat sink made of NePCM. According to their experimental results, dispersing 3% of the nanoparticles in pure PCM increases the viscosity and thermal conductivity of NePCM by 100% and 150%, respectively. They discovered that adding 0.5% of the nanoparticles to pure PCM could achieve a temperature reduction of 15 • C for the heat sink. Furthermore, their findings confirmed that when nano-additives were added in concentrations greater than 0.5%, the thermal efficiency of heat sinks decreased significantly. In another experiment, Li et al. [56] evaluated PCM's heat transfer and energy storage properties containing graphene nanoplatelets (GNP) while charging it in a rectangular cavity. They discovered that increasing the GNP concentration decreases both heat transfer rate and thermal storage simultaneously. They claimed that the increase in thermal conductivity caused by the addition of GNP was insufficient to compensate for PCM's increase in viscosity, which resulted in a significant decrease in natural convection. They stated that prior research had grossly underestimated the detrimental effects of increased viscosity caused by the addition of nanomaterials. Praveen et al. [57] investigated the heat transfer performance of graphene nano-platelets laden microencapsulated PCM (ME/GNP PCM) in a finned thermal energy storage-based heat sink. They added GNP of 0.5, 1, and 3 wt% to enhance the heat transfer between the encapsulated PCM. Their results showed that the GNP in the encapsulated PCM reduced the recovery time of the heat sink. Ramakrishnan et al. [58] examined the impacts of different carbon additives on the thermal storage performance of form-stable PCM integrated cementitious composites. They utilized various high conductive carbon-based additives such as graphite (G), carbon nanotubes (CNT), and graphene nanoplatelets. Their outcomes confirmed that although G and GNP additives showed a high thermal conductivity increment, heat transfer performance tests proved that GNP leads to the highest performance enhancement, and graphite leads to the least. Singh et al. [59] analyzed the thermal enhancement of a binary eutectic PCM, laden with different concentrations of COOH-functionalized graphene nanoplatelets (f-GNP) for a multi-effect solar cooling thermal storage system. They investigated the effects of doping f-GNP on the conjugate heat transfer inside the PCM and fluid flow of HTF in a vertical shell and tube type storage system, suitable for the double effect solar cooling system. Their findings illustrated that the f-GNP dispersions accelerate the heat storage process, with a maximum decrease of 17.3% in the total melt duration. Mahdi et al. [60] investigated the melting augmentation in a triplex-tube latent heat energy storage system using a nanoparticle-metal foam combination. They claimed that dispersing nanoparticles in the presence of metal foams results in melting time savings of up to 90% depending on the foam structure and volumetric nanoparticle concentration. They observed that the melting time decreases as the porosity decreases and/or volume fraction increases. They recommended using high-porosity metal foam with low volume-fraction nanoparticles to promote the positive contribution of natural convection during the melting process. Senobar et al. [61] established an experimental setup to study the charging/discharging performance of PCM-nanoparticle-metal foam composites. They showed that adding nanoparticles to PCM enhanced the heat transfer rate by up to 24%. Moreover, PCM-metal foam revealed a heat transfer rate that was up to 26% faster when compared with pure PCM. They found that by using metal foam and nanoparticles simultaneously, the heat transfer rate increased significantly (up to 65% compared to pure PCM).

Shell and Tube LHTES Systems
As previously stated, solar energy is a time-dependent energy source, which means there is no match between energy supply and demand, and this is a significant disadvantage. In this case, combining solar energy systems and LHTES units may be an effective solution. As a result, TES systems can be used to rationalize the use of solar thermal energy in a variety of applications, including desalination [62], distillation [63], drying [64], industrial process heat [65], hot water [66], and air conditioning [67]. In most solar energy applications, a shell-and-tube TES system is used in which the heat transfer fluid (HTF) flows inside the tube, and the PCM is installed on the shell side. The PCM is subjected to the melting/solidification process via HTF temperature variations. Adopting techniques that accelerate the charging/discharging process and increase the amount of stored energy in shell-and-tube LHTES systems can be highly beneficial. Numerous researchers have attempted to improve the efficiency of shell-and-tube LHTES systems, and several of these efforts are summarized below.
Masoumi and Khoshkhoo [68] studied the impacts of simultaneously installing fins and dispersing nanoparticles on the duration of PCM charging in a horizontal shelland-tube TES unit. Their findings indicate that adding 0.39 wt% titanium dioxide (TiO 2 ) nanoparticles to stearic acid (SA) as PCM serves to reduce charging duration by approximately 4%, while fins reduce charging duration by 68%. They stated that adding nanoparticles reduces the effectiveness of natural convection as the primary heat transfer mechanism in PCM charging. Shi et al. [69] simulated the charging/discharging process of PCM in a shell and multi-tube LHTES system under the influence of magnetic fields as an active enhancement technique in a numerical study. They discovered that applying a magnetic field enhanced free convection in molten PCM. Additionally, increasing the magnetic field's intensity results in a decrease in charging/discharging time. They claimed that charging and discharging times were reduced by 80.02% and 53.19%, respectively, in the presence of a magnetic field. Karami and Kamkari [70] used perforated fins in their empirical work to attempt to improve the natural convection in vertical shell and tube LHTES systems. They used lauric acid and water as PCM and HTF, respectively. They contended that perforating the fins would improve natural convection within the PCM shell. They discovered that the charging time was 7% shorter for the perforated finned shell and tube unit than for the original unit with solid fins. Mahdavi et al. [71] investigated the charging and discharging efficiency of a vertical shell and tube LHTES unit through the use of heat pipes and nano-additives. The heat pipes were used to connect the PCM tank to the HTF. They dispersed the nanoparticles of four metals and metal oxides in PCM, including Al 2 O 3 , Ag, Cu, and CuO. They demonstrated that the addition of nanoparticles significantly reduced the charging/discharging time and energy stored in the system. They discovered that a particular type of nanoparticle effectiveness depends on the system's number of embedded heat pipes. Mahood et al. [72] examined the effect of fin arrangement on the thermal efficiency of a finned shell and tube LHTES system using numerical simulations. They discovered that increasing the fin length to 0.8 of the shell radius significantly reduced the charging duration by 50%. Additionally, they confirmed that changing the angle of the fins from 72 • to 15 • increases the melting time by 40%.

Research Contribution
The total charging/discharging time and the total stored energy are critical design parameters for LHTES units. The various heat transfer augmentation methods described in the literature have a major impact on these two factors. For example, while the use of fins accelerates the phase change process in LHTES systems, it also reduces the total stored energy, which is a disadvantage [72][73][74]. Numerous previous studies have recommended the utilization of nano-additives as an effective technique for increasing the efficiency of LHTES units. However, it should be noted that the dispersion of nano-additives can significantly alter the properties of PCMs, such as density, specific heat, latent heat, viscosity, and thermal conductivity. These substantial differences can result in unpredictable dynamic and thermal behaviors for various PCMs. Thus, when nano-additives are used in LHTES systems, their thermophysical properties (thermal conductivity, density, and specific heat capacity) and concentrations are critical. The importance of this issue, which has been overlooked in previous research, motivated the current study to investigate the effects of the concentration and thermophysical properties of nano-additives on phase change duration and total stored energy, which are used as basic parameters in the design of LHTES units. To this end, an optimization-based procedure was required, which begins with a thorough examination of nano-additives' effects on the efficiency of the shell and tube LHTES unit and ends with the principled adoption of optimal solutions.
The present research employed a hybrid approach based on the design of experiments (DOE), computational fluid dynamics (CFD), artificial neural networks (ANNs), multiobjective optimization (MOO), and multi-criteria decision making (MCDM) to optimize the concentrations and properties of nano-additives dispersed in a shell and tube LHTES unit containing paraffin wax as PCM.
Five steps comprised the proposed hybrid procedure. The first step is to determine the total number of numerical runs required for the optimal search of the design space using a DOE technique and the range of optimization variables. The second step simulates the charging of NePCM in a two-dimensional shell and tube LHTES unit using CFD for samples provided by DOE in order to determine the effect of optimization variables (thermal conductivity, density, specific heat, and nano-additive concentrations) on the objective functions (total charging time and stored energy). The third step, which utilizes ANNs, presents polynomial models for objective functions in terms of optimization variables. MOO is applied to the obtained mathematical models in the fourth step to generate a set of optimal Pareto points. Finally, selected optimal points with various characteristics are presented using various MCDM methods in the fifth step. Figure 1 depicts the proposed hybrid procedure's flowchart in greater detail. The approach presented in this paper as an effective method for selecting NePCM for shell and tube LHTES systems can pave the way for multiple researchers, particularly those working on solar energy applications.   Figure 2 illustrates an overview of a vertically oriented shell and tube LHTES unit in a solar thermal storage system. As can be seen, the LHTES charging process begins as the solar collector's heated HTF flows upward through the inner tube. Heat transfer from the HTF to the shell increases the temperature of the PCM, which continues until the PCM is completely melted. To use the stored energy at the appropriate time, the HTF descends through the inner tube where the thermal energy is transferred to the HTF via a solidification process that increases the HTF's temperature. In our research, four basic characteristics of nano-additives were considered as optimization variables: concentration and thermophysical properties such as specific heat, density, and thermal conductivity. The volume fraction of nano-additives (the volumetric concentration of the nanoparticles in the PCM) is a critical parameter affecting engineering systems' dynamic and thermal behavior. Most previous research has demonstrated that increasing the volume fraction of nano-additives (ϕ) above a certain point reduces the system's efficiency [75,76]. On the other hand, the models used to predict the thermophysical properties of NePCM will be inaccurate for ϕ > 0.05 [74]. As a result, the maximum volume fraction of nano-additives was set to 0.05. Thermal conductivity, density, and specific heat properties of nano-additives that are widely utilized in the industry can be used to determine the study range for thermal conductivity, density, and specific heat. The thermophysical properties of common nano-additives used in engineering systems are summarized in Table 1 [77][78][79][80][81][82][83]. Therefore, as per Table 1 and the volume fraction of nano-additives considered, the variation range of optimization variables is presented in Table 2. Table 1. Thermophysical properties of widely used nano-additives [77][78][79][80][81][82][83].

Optimization Variables Unit From To
Thermal conductivity (k) By employing a proper design of experiments method, the design space can be analyzed optimally. To this end, the well-known D-optimal method [84,85] was used to minimize calculation costs and extract as much information from the variables' space as possible. Through the D-optimal design technique, 45 DOE points were generated and distributed within the design space. These DOE points are used to populate the CFD model with input parameters.

Assumptions and Boundary Conditions
The following are the assumptions that govern the current problem:

•
It is assumed that a fully developed turbulent HTF flow enters the LHTES system.
The following equation provides an estimate of the entrance length required to create a fully developed turbulence flow [86]: The inner tube has a wall thickness of 2 mm and is made of aluminum.

•
The thermophysical properties of the PCM are temperature-dependent.

•
The thermophysical properties of the HTF are considered constant.

•
It is assumed that the nano-additives are dispersed homogeneously in the PCM.

•
The liquid NePCM flows in the LHTES unit are assumed to be Newtonian, steady-state, incompressible, and laminar flows.
The boundary and initial conditions can be expressed as follows: • The Reynolds number and inlet temperature of HTF are set to be 4500 and 363.15 K, respectively.

•
The outer surfaces of the LHTES system are considered perfectly insulated (adiabatic condition). • A zero gradient boundary condition (outflow) is considered for the HTF outlet. • A 2D axisymmetric modeling is adopted for the simulation of the LHTES system (see Figure 2).

•
No-slip condition is assumed at the boundary of the walls.

•
The initial temperature of the NePCM is set to be 310.15 K (the NePCM is subcooled).

Thermophysical Properties of NePCM
Paraffin wax was used as a common organic PCM to store thermal energy in this research as it is one of the most widely used PCMs in mid-temperature TES applications, such as building air conditioning, solar water distillation, solar collector systems, and electronic thermal management [87][88][89][90]. Paraffin wax has a wide range of applications due to its desirable properties. This PCM is odorless, chemically inert, non-toxic, and inexpensive [91]. During the phase change process, it has a low volume change and a low vapor content in melt [92]. Additionally, paraffin wax is non-reactive to most encapsulation materials, is compatible with almost all materials, has a high latent heat, and exhibits chemical stability over numerous charging/discharging cycles [93]. However, its primary disadvantage is its low thermal conductivity, which can be compensated for through enhancement techniques [94]. Furthermore, water is used as an HTF in conventional solar collector systems. Table 3 presents the thermophysical properties of paraffin wax and water [95,96]. Table 3. Thermophysical properties of PCM and HTF [95,96].

Properties Unit Paraffin Wax (PCM) Water (HTF)
Specific heat, In most previous studies, the thermophysical properties of nano-enhanced PCMs have been estimated using analytical models and experimental correlations. Fan et al. [97] tested several conventional models for predicting the properties of NePCM in their experimental study on the performance of TES systems. They observed that when the concentration of nano-additives dispersed in PCM is less than five vol%, conventional models and experimental correlations report an acceptable error compared to the experimental values. Thus, using simple models based on mixing theory [98] to calculate the properties of NePCM, such as density, latent heat, thermal expansion, and specific heat, can be considered reliable over a certain range of nano-additive concentrations. The reduced form of the mixing model is employed as follows: Vajjha et al. [99] proposed a model for estimating the viscosity of NePCM that is reasonably accurate: The thermal conductivity of NePCM containing uniformly sized, spherical nanoparticles is calculated using the well-known Maxwell model [100]: Several previous numerical studies in various fields have used the expressed models to calculate the thermophysical properties of NePCM [101][102][103][104][105][106].

Governing Equations and Numerical Procedure
Based on assumptions, transient continuity, momentum, and energy equations are introduced as follows [107]: where ρ denotes the density, H denotes the enthalpy, k is the thermal conductivity, µ denotes the dynamic viscosity of liquid PCM, P denotes the pressure, and T re f represents the reference temperature. The solid-liquid phase change process is simulated and tracked using an enthalpyporosity method. Voller and Prakash [108] described this technique as one that does not require explicitly tracking the phase's interface. The enthalpy-porosity method has been used in recent research to improve the accuracy of phase change simulations [109][110][111][112].
Source term → S in Equation (9) represents the impact of the enthalpy-porosity technique and is defined by the following equation: where A mushy denotes a mushy zone constant and estimates the amplitude of damping [113].
Although the mushy zone constant varies between 10 4 and 10 7 , the traditional value of 10 5 has been used in most previous studies using conventional PCMs [74,107]. Additionally, ε by a small constant value of 0.001 is introduced to avoid division by zero. Moreover, λ is a critical parameter in phase change analysis because it indicates the liquid fraction of PCM; λ possesses a value between 0 and 1 that is determined by the PCM temperature: where T l and T s are the liquidus and solidus temperatures, respectively. Natural convection takes precedence over thermal conduction once the melting process begins and the liquid level is raised and is recognized as the main heat transfer mechanism in the TES system [74,107]. Consequently, the well-known Boussinesq approximation is employed to model the influence of natural convection (buoyancy force) as a melting accelerator mechanism. The Boussinesq model calculates the buoyancy term by considering density as a function of temperature [114]: However, in calculating other density-dependent terms, the density value is considered constant.
The enthalpy term H in Equation (10) is calculated as follows: where ∆H and h denote the latent and sensible enthalpy, respectively, and can be obtained as: where h re f is enthalpy at the reference temperature (T re f ). The realizable k − ε turbulence model with a standard wall function was selected to simulate turbulent HTF flow in the inner tube.
The present two-dimensional axisymmetric model was constructed using an ANSYS Fluent CFD solver through the finite volume method. For pressure-velocity coupling, the well-known SIMPLE algorithm was utilized. The QUICK scheme is used to discretize both the energy and momentum equations. Additionally, a second-order implicit transient formulation is used for increased precision, and the PRESTO pressure discretization scheme is implemented [115]. Under-relaxation parameters for energy, liquid fraction, momentum, and pressure were set to 1.0, 0.9, 0.7, and 0.3 to reduce fluctuations and increase stability during the solving process.
Moreover, to ensure a convergence of results, the convergence criteria for the momentum and continuity equations were adjusted to 10 −4 , while the convergence criteria for the energy equations were adjusted to 10 −8 . Several runs were conducted as a trial-and-error strategy to determine the most appropriate time step size for ensuring transient solution stability; finally, the optimum time step size of 0.1 s was obtained. The number of iterations per time step was set to 100 to ensure that each time step was convergent.

Grid Sensitivity Study and CFD Validation
Using ICEM software, the computational domain for the shell and tube LHTES system was drawn and meshed. The structured meshing technique was applied to generate a quality grid. Furthermore, greater emphasis was placed on the quality of the grids in the NePCM zone and near the wall boundaries to perfectly capture gradients. Five grid designs with 9505, 14,617, 23,213, 36,118, and 50,028 nodes were analyzed to determine the CFD results' independence from the grid size. This was accomplished by comparing the transient liquid fraction and the average temperature of pure PCM across all grids, as shown in Figure 3. As can be seen, the maximum change in the outcomes was minimal once the grid size was increased to 36118. The grid sensitivity procedure was repeated for NePCM cases, and the results were comparable. As a result, the grid structure with 36118 nodes was deemed adequate for examining all cases in the current work. Figure 4 displays the solution domain and selected grid structure in the present research.
The precision and reliability of the developed computational model were evaluated using the experimental and numerical results of Seddegh et al. [116] and the experimental data tabulated by Akgun et al. [117]. Both studies used HTF to flow through the inner tube and paraffin wax to fill the shell space, identical to the current TES system. Seddegh et al. [116] used visualization to investigate the solid-liquid interface during the PCM charging process in a vertical cylindrical shell-and-tube LHTES system. They used numerical simulations to elucidate and expand on their experimental results. Figure 5a compares the present numerical results to the experimental and numerical data of Seddegh et al. [116] for the transient volumeaverage temperature of PCM.   On the other hand, in the study by Akgun et al. [117], the local temperature was determined using a total of 24 thermocouples located within the PCM shell. The transient temperature of selected points (T13 and T61) was used to compare the present results to their experimental data. T13 was located at a radius of 47 mm and a height of 0 mm, while T61 was located at a radius of 28 mm and a height of 421 mm, with radial coordinates taken from the inner pipe's center. Figure 5b illustrates the details of the findings' comparisons.
The following statistical criteria, namely mean absolute percentage error (MAPE), were used to quantify the comparison between the current numerical outcomes and the findings of [116,117]: where Y i and X i denote previous and present findings, respectively. Table 4 compares the current results to previous studies that used MAPE criteria. As can be seen, MAPE was less than 5% in all four comparisons, confirming the validity of the current numerical model for examining the PCM melting process in shell and tube LHTES systems.

Objective Functions
Two critical factors in the scope of TES that have dominated researchers' efforts in designing LHTES systems are the duration of energy storage and the amount of stored energy. Several authors have concentrated on charging/discharging times [118][119][120], while others have concentrated on total stored energy [121,122]. Meanwhile, a few authors have concentrated on optimizing these two factors concurrently [112].
The present study simultaneously evaluated both total stored energy and melting duration as objective functions during the optimization process. In fact, changes in the volume fraction and properties of nano-additives as optimization variables have a complex effect on objective functions. These effects were investigated using CFD simulations, and, finally, the optimal conditions could be presented through a robust optimization process. The melting time and total stored energy in the LHTES system were calculated using the following plan: where λ vol−avg denotes the volume-average liquid fraction and t m is melting time. In the present study, the total stored energy (E t ) refers to the total stored energy at t m .

GMDH-Type Neural Network Modeling
Artificial neural networks (ANNs) are a powerful tool for modeling and forecasting various highly complex phenomena [123][124][125][126][127][128][129][130][131][132][133]. They are inspired by the biological neural networks found in the human brain. Even though various ANN types are capable of modeling a wide variety of systems, their accuracy and efficiency are dependent on the quantity and quality of available data. Ivakhnenko [134] developed a group method of data handling (GMDH) polynomial neural networks with a self-organizing algorithm and a structure similar to feed forward ANNs to make modeling accuracy less dependent on the quantity of data. In the GMDH method, the modeling process begins with the input layer, which comprises the input data neurons. Quadratic polynomials connect pairs of neurons from various layers to form neurons in the middle layers, and a single neuron is presented as a predicted model of the output data in the output layer. One of the most remarkable properties of the GMDH algorithm is its self-organizing process in which various configurations of pairs of neurons are examined and only those that contribute to the output model's efficiency are retained in the middle layers.

Multi-Objective Optimization
The optimal conditions for the optimization variables should result in a shorter melting time and an increase in the amount of energy stored in the LHTES system. Due to the complexity of the objective functions' behavior in the presence of varying nano-additive properties, it is necessary to specify conditions for the variables used to optimize the objective functions. For this purpose, the Pareto-optimal procedure is applied. A set of solutions (Pareto front) is presented as the optimal points in this procedure. All optimal Pareto points are equally desirable, and the designer's strategy determines their selection. The Pareto front for the current optimization problem was calculated using the well-known TOPSIS algorithm [142]. Numerous designers prioritize this algorithm when solving multi-objective optimization problems, and its efficiency has been compared to that of other methods in various fields [143][144][145][146][147]. The NSGA-II optimization process is summarized in Figure 8. The two-objective optimization problem in this work can be formulated as follows: Figure 9 illustrates a comparison of Pareto front values with CFD results. As can be seen, when ANN and MOO were combined, the optimal set of solutions could be found in comparison to the original CFD data. (50) Figure 9 illustrates a comparison of Pareto front values with CFD results. As can be seen, when ANN and MOO were combined, the optimal set of solutions could be found in comparison to the original CFD data.  Figures 10 and 11 show the optimal variation of melting time and stored energy in terms of optimization variables for all optimal points. As depicted in Figures 10a and 11a, for all Pareto points in the range of 0.5-4.9%, ϕ was uniformly distributed. In other words, the simultaneous optimization of t m and E t is not restricted to specific values and can be achieved by adding any volume fraction of nano-additives. However, as shown in Figures 10a and 11a, phase change time and stored energy increased as the volume fraction of nano-additives increased. This means that if the melting duration is more critical, a lower concentration of nano-additives should be used; however, if the stored energy is more critical, the volume fraction of nano-additives should be increased. The volume fraction trend described above was not observed for the thermophysical properties. According to Figures 10a and 11a, approximately 97% of optimal points lay within the range 180 ≤ k NP ≤ 480. By increasing the thermal conductivity of nano-additives, the probability of being one of the optimal solutions increased to the point where more than half of the total Pareto points were in the range 315 ≤ k NP ≤ 480. In general, increasing the thermal conductivity of nano-additives can reduce phase change time and stored energy in the TES system. additives should not be less than a particular value; this reduces the system's performance. Around 86% of the optimal points lay between 9950 < < 9990 , indicating that the dispersion of high-density nano-additives has a beneficial effect on the paraffin wax melting process. The trend described above was more pronounced for . This means that, for approximately 95% of the optimal points , the values exceeded 795. Consequently, TES systems can significantly benefit from the addition of nanoadditives with high specific heat.
(a) (b) Figure 10. Optimal variation of melting time ( ) with respect to (a) volume fraction and thermal conductivity, (b) specific heat and density. The results of the Pareto front show that the dispersion of nanoparticles (especially with a high concentration) in paraffin wax does not have a favorable effect on the phase change time. On the one hand, raising the volume fraction of nanoparticles enhances the thermal conductivity of NePCM and, on the other hand, increases the viscosity of NePCM. Growing the viscosity of NePCM diminishes the impact of natural convection as the primary mechanism of heat transfer in the PCM melting process. In fact, increasing the nanoparticles concentration has a two-way impact. On the one hand, it raises the conduction heat transfer; on the other hand, it causes a decrease in convection heat transfer, which subsequently increases the phase change time. To clarify this point, it should be noted that using pure PCM takes 2245 s for complete melting. This is while in the presence of nanoparticles; therefore more time is needed for the charging process. In addition, increasing the concentration of nanoparticles affects the amount of total stored energy. As mentioned, the total thermal energy stored in LHTES systems includes latent and sensible heat energy. According to Equation (4), by increasing the volume fraction of nanoparticles, the amount of latent energy decreases, however, on the other hand, due to the high density of nanoparticles, the amount of sensible energy improves. By raising the volume fraction of nanoparticles, the increment in the sensible energy overcomes the reduction in latent heat energy, and these interactions increase the total stored energy in the TES system. As confirmation of this viewpoint, it can be stated that the total amount  Figures 10b and 11b show a similar trend for the density and specific heat of nanoadditives, indicating that the congestion of the optimal points within a specific range of ρ NP and C p NP is high and implying the optimal range of these properties for selecting the appropriate nano-additives. As can be seen, all optimal solutions had ρ NP values between 9740 and 9990. No optimal point exists when ρ NP < 9740. As a result, it is essential to note that, when studying TES systems containing paraffin wax, the density of dispersed nano-additives should not be less than a particular value; this reduces the system's performance. Around 86% of the optimal points lay between 9950 < ρ NP < 9990, indicating that the dispersion of high-density nano-additives has a beneficial effect on the paraffin wax melting process. The trend described above was more pronounced for C p NP . This means that, for approximately 95% of the optimal points C p NP , the values exceeded 795. Consequently, TES systems can significantly benefit from the addition of nano-additives with high specific heat.
The results of the Pareto front show that the dispersion of nanoparticles (especially with a high concentration) in paraffin wax does not have a favorable effect on the phase change time. On the one hand, raising the volume fraction of nanoparticles enhances the thermal conductivity of NePCM and, on the other hand, increases the viscosity of NePCM. Growing the viscosity of NePCM diminishes the impact of natural convection as the primary mechanism of heat transfer in the PCM melting process. In fact, increasing the nanoparticles concentration has a two-way impact. On the one hand, it raises the conduction heat transfer; on the other hand, it causes a decrease in convection heat transfer, which subsequently increases the phase change time. To clarify this point, it should be noted that using pure PCM takes 2245 s for complete melting. This is while in the presence of nanoparticles; therefore more time is needed for the charging process. In addition, increasing the concentration of nanoparticles affects the amount of total stored energy. As mentioned, the total thermal energy stored in LHTES systems includes latent and sensible heat energy. According to Equation (4), by increasing the volume fraction of nanoparticles, the amount of latent energy decreases, however, on the other hand, due to the high density of nanoparticles, the amount of sensible energy improves. By raising the volume fraction of nanoparticles, the increment in the sensible energy overcomes the reduction in latent heat energy, and these interactions increase the total stored energy in the TES system. As confirmation of this viewpoint, it can be stated that the total amount of stored energy in the case of pure PCM is 251111 Joule. However, in all optimal Pareto points, a higher amount of total thermal stored energy is observed.

Multi-Criteria Decision Making
The NSGA-II algorithm was used to determine Pareto optimal points. Pareto points are a collection of all optimal solutions that do not overlap but remain superior to all other solutions in the problem space. This means that, in the current problem with two objective functions, changing the optimal solution does not improve both objective functions simultaneously but instead alters their importance to one another, which is where the designer can use MCDM techniques to select the desired optimal points from numerous Pareto front solutions. The well-known TOPSIS [148] and VIKOR [149] methods are employed for this purpose. These two methods work based on the exact mechanism of aggregation function (closeness to the ideal point). However, the criteria are normalized by two different techniques. The TOPSIS method uses Euclidean normalization while the VIKOR method employs linear normalization [149]. Numerous previous studies have utilized these two techniques and demonstrated their efficacy and effectiveness [32,[150][151][152]. The following describes the process of decision making using the TOPSIS and VIKOR methods: A worst = max t ij i = 1, 2, . . . , m j ∈ J − , min t ij i = 1, 2, . . . , m j ∈ J + ≡ t wj j = 1, 2, . . . , n where 5 Computing the alternative Euclidean distance for the best and worst alternatives: Ranking the alternatives from 0 to 1 as follows: The alternative with the highest S i,worst value reaches the best rank in the TOPSIS method.

VIKOR
The alternatives' ranking process by the VIKOR technique can be described as follows [112]: 1.
Forming the decision-making matrix (x ij ) m×n 2.
Obtaining the best and worst alternatives for all criteria.
Computing the regret measure ( R i ) and utility measure ( S i ): 2.
Computing the VIKOR index (Q i ): The value of v is set to 0.5 to create a moderate strategy between regret and utility measures. In the VIKOR method, the alternative with the lowest Q i value achieves the best rank.

Selection of Optimal Points
If w t and w E denote the importance (weight) of the t m and E t , respectively, there are three strategies for selecting the desired points: Unbalanced weights (w t > w E or w t < w E ). Table 6 summarizes the optimal points determined using the strategies above for the TOPSIS and VIKOR methods. As shown, the TOPSIS method provides points A-C for various weights. Points D-F are suggested by the VIKOR method. Additionally, G and H points demonstrate single-criteria decision making for t m and E t , respectively. Moreover, for further clarity, Figure 12 illustrates the location of the selected points on the Pareto front. As can be seen, the proposed points for the same weights by the TOPSIS and VIKOR techniques are not identical. For example, when w t = w E = 0.5, the TOPSIS method recommends point A, whereas the VIKOR method recommends point B, which is a significant difference. This trend is stronger when w t = 0.25 and w E = 0.75. The question now is which method is the superior option. The deviation index (DI) can be used to address this question [153]. This index, calculated as follows [154], assigns a value to alternatives based on their position relative to the ideal and non-ideal solutions (see Figure 12).
where t m n and E t n represent Euclidian, non-dimensioned t m and E t , respectively. According to Equations (67)- (69), the solution closer to the ideal point and farther away from the non-ideal point has a lower DI value and is, therefore, more appropriate. The DI values for the TOPSIS and VIKOR methods are shown in Table 7. As can be seen, when the objective weights are w t = 0.75 and w E = 0.25, the VIKOR method provides a superior solution; however, for other objective weights, the TOPSIS method's proposed solution takes precedence. Mathematics 2021, 9,    According to Equations (67)- (69), the solution closer to the ideal point and farther away from the non-ideal point has a lower DI value and is, therefore, more appropriate. The DI values for the TOPSIS and VIKOR methods are shown in Table 7. As can be seen, when the objective weights are = 0.75 and = 0.25, the VIKOR method provides a superior solution; however, for other objective weights, the TOPSIS method's proposed solution takes precedence. In order to ensure the accuracy of the combined optimization process in this study, selected points A-H were numerically re-simulated. Table 8 shows that the maximum error of optimization outcomes with CFD results were less than 0.5%, indicating that the combined optimization process has acceptable precision.  In order to ensure the accuracy of the combined optimization process in this study, selected points A-H were numerically re-simulated. Table 8 shows that the maximum error of optimization outcomes with CFD results were less than 0.5%, indicating that the combined optimization process has acceptable precision. It should be noted that the weight of the objective functions to select the suitable solution from the Pareto points (MCDM process) can depend on various issues such as environmental conditions, designer strategy, time and frequency of system operation, and other constraints. Consider a solar energy storage system. If this system is limited in time and number of phase change cycles, an enormous weight should be assigned to t m . On the other hand, if the time and number of cycles are not restricted, the amount of stored energy should be prioritized, i.e., more weight should be assigned to E t .

Conclusions
One of the primary concerns of heat transfer engineers has been the application of phase change materials containing nano-additives with optimal properties that increase stored energy and decrease phase change time in TES systems. The present paper employs a combined procedure based on the design of experiments (DOE), computational fluid dynamics (CFD), artificial neural networks (ANN), multi-objective optimization (MOO), and multi-criteria decision making (MCDM) to optimize the properties of nano-additives dispersed in a shell and tube LHTES system containing paraffin wax as a phase change material (PCM).
Five steps comprise the proposed hybrid procedure. The first step is to determine the total number of numerical runs required to conduct an optimal search of the design space using a DOE technique. The second step employs the CFD method to analyze the melting process in the shell and tube LHTES system. In the third step, polynomial models for objective functions are presented in terms of optimization variables using the GMDH method. MOO is used in the fourth step to generate a set of optimal Pareto points. Finally, in the fifth step, various MCDM methods are used to introduce optimal points with various characteristics.
The following summarizes the optimal variation of objective functions: • ϕ is uniformly distributed for all Pareto optimal points in the range of 0.5-4.9%. • By increasing ϕ in optimal solution, the melting time and stored energy are enhanced. • Approximately 97% of optimal points lie within the range 180 ≤ k NP ≤ 480. Thus, the addition of nanoparticles with high thermal conductivity can increase the efficiency of TES systems. • Around 86% of the optimal points lie between 9950 < ρ NP < 9990, indicating that the dispersion of high-density nano-additives has a beneficial effect on the paraffin wax melting process. • Approximately 95% of optimal points have C p NP values greater than 795. Consequently, TES systems can benefit significantly from the addition of nano-additives with high specific heat.
Proper validation should be performed at each stage of the hybrid technique to ensure the accuracy of the final solutions. In fact, ensuring the high accuracy of CFD, ANN, MOO, and MCDM methods is exceptionally significant and can guarantee the reliability of optimal solutions. However, the sequence of the various steps is also very effective in the output results, in addition to accuracy. This means that removing DOE from the hybrid technique causes the problem space to be inadequately explored, resulting in a decrease in the efficiency of GMDH output models. Furthermore, the low precision of the CFD results and improper settings of GMDH and NSGA-II schemes cause inaccuracies in the final output of the problem. In future work, by applying the present technique, the geometric characteristics of a shell-and-tube LHTES unit will be optimized.