DOI: 10.3303/CET2188050 Paper Received: 3 May 2021; Revised: 24 September 2021; Accepted: 2 October 2021 Please cite this article as: Zálešák M., Charvát P., Klimeš L., 2021, Robustness and Accuracy of the Particle Swarm Optimisation Method in the Solution of Inverse Heat Transfer Problems with Phase Change, Chemical Engineering Transactions, 88, 301-306 DOI:10.3303/CET2188050 CHEMICAL ENGINEERING TRANSACTIONS VOL. 88, 2021 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Petar S. Varbanov, Yee Van Fan, Jiří J. Klemeš Copyright © 2021, AIDIC Servizi S.r.l. ISBN 978-88-95608-86-0; ISSN 2283-9216 Robustness and Accuracy of the Particle Swarm Optimisation Method in the Solution of Inverse Heat Transfer Problems with Phase Change Martin Zálešáka,*, Pavel Charváta, Lubomír Klimešb a Department of Thermodynamics and Environmental Engineering, Brno University of Technology, Technická 2896/2, 61669 Brno, Czech Republic b Sustainable Process Integration Laboratory—SPIL, NETME Centre, Brno University of Technology, Technická 2896/2, 61669 Brno, Czech Republic 162098@vutbr.cz Some engineering heat transfer problems cannot be solved directly as certain parameters influencing heat transfer are unknown. This is where the solution of inverse heat transfer problems might be utilized. The paper explores application of the Particle Swarm Optimisation (PSO) method to heat transfer problems with phase change. In order to assess the robustness and accuracy of the method, the study was conducted on simulation results. The explored case is an air-PCM heat exchanger. In the first step, the simulations with various values of parameters influencing heat transfer were performed. All cases were simulated as transient with time- dependent boundary conditions. In the second step, some of these parameters (phase change temperature, density, thermal conductivity, etc.) were assumed to be unknown and they were obtained from the simulation results with the use of the PSO method. The cost function was defined in the form of a root mean square error between the simulation results (the outlet air temperature from the heat exchanger) and the results of the optimized scenario. The number of unknown parameters varied from one to seven (in case of the description of the relationship of effective heat capacity on temperature, which was parameterized in the form of an asymmetric Gaussian function). The results have shown that the PSO is a robust and accurate tool with only up to 2 % difference in the enthalpy of fusion and no significant discrepancies during the higher dimension optimisation. 1. Introduction The excessive usage of traditional fossil fuels in order to cover the majority of daily energy demands may lead to many emerging environmental concerns such as the energy shortage crisis, global warming and environmental pollution (Liu et al., 2021). The development and utilisation of renewable energy as a way to at least partially replace fossil fuels becomes more and more significant (Gossard et al., 2015). In recent years, there has been an increasing usage of Phase Change Materials (PCMs) as the heat storage medium in the Thermal Energy Storage (TES) systems. The main advantage of PCMs are a higher energy density (due to the latent heat) and improved overall efficiency of TES. PCMs are as well very useful as a TES medium in cases where there is not a clear overlap between the energy supply and energy demand e.g. in case of the solar radiation (Grabo et al., 2019). In order to assess the potential of PCMs in the industrial applications, the development of the corresponding numerical model is needed. Such a model always requires input parameters such as the material properties, boundary conditions and design parameters. The overall accuracy of a numerical model highly depends on the thermal material properties, which have to be measured with a suitable experimental method. A wide range of experimental methods is used for these purposes e.g. differential scanning calorimetry (DSC), differential thermal analysis (DTA) or the temperature history (T-history) method (Cascone and Perino, 2015). All mentioned techniques are based on the monitoring of the thermal behaviour of a sample of the studied material and comparing it with the material with known thermophysical properties. The main drawback of the DSC, which is the most common experimental technique for characterization of PCMs, is the sample size (usually up to 10 mg 301 or less), which may lead to inaccuracies since the most of practical applications are rather large scale (Gossard et al., 2015). As an alternative approach, the thermophysical properties can be determined as a solution of the inverse heat transfer problem. A study dealing with the identification of the enthalpy of phase change using the inverse problem and calorimetry experiments was conducted by (Franquet et al., 2012). The authors have developed the numerical model of a cylindrical-shape sample of PCM to be used for inverse identification. The objective function was constructed in terms of the difference between the numerically simulated and experimentally measured heat fluxes. As an optimisation solver, the Nelded-Mead simplex optimisation method was used. The proposed method was validated using several different samples of material. The experimentally measured enthalpy curves were compared to the enthalpy curves obtained as a solution of the inverse problem. The authors have concluded that the proposed method estimated the behaviour well and it has shown to be very accurate. Cascone and Perino (2015) have presented the method for estimation of the specific heat-temperature dependency of PCMs employing the inverse problems. The experimental measurements of the sample surface temperatures were made under boundary conditions and the inverse problem was formulated in order to minimize the difference between the simulated and measured surface heat fluxes. The DSC measurements were used for compassion of the specific heat curves. The results have shown that the peak value of the specific heat curve and the specific heat in the liquid state were comparable with the melting curve obtained by DSC, on the other hand the estimation of phase change temperature was congruent with the DSC solidification curve. The present study is focused on the robustness and accuracy of the PSO method used for the inverse identification of the material properties of a paraffin-based PCM. The authors are following up on the previous study (Zálešák et al., 2020), where the PSO was adapted during the design optimisation of a solar air collector. 2. Methods The air-PCM heat exchanger (air-PCM HX) was considered in this study (see Figure 1a) as a thermal system for inverse identification of PCM properties. The air-PCM HX consisted of 100 Compact Storage Modules (CSMs); aluminium plates filled with paraffin-based PCM. In practice, incorporating CSMs is a common and practical way of PCM integration in TES applications, since the desired overall heat storage capacity of the TES system can be easily achieved by a suitable arrangement of the appropriate number of CSMs. The dimensions of the CSMs considered in this study were 450 mm x 300 mm x 10 mm and they were arranged in 5 rows of 20 CSMs. Figure 1: (a) The air-PCM HX and (b) the corresponding schematics of the unit 2.1 Numerical model The proposed numerical model of the air-PCM HX was constructed as a quasi-2D heat transfer problem. The numerical model consists of a series of 1-D heat transfer problems (heat conduction through a PCM layer in the direction of its thickness) and the governing equation could be formulated as 𝜌ceff 𝜕T 𝜕𝜏 = 𝜕 𝜕𝑥 (𝑘 𝜕T 𝜕𝑥 ) (1) where T is temperature, 𝜌 is density, 𝜏 is time, 𝑘 is thermal conductivity, ceff stands for effective heat capacity and x is the spatial coordinate in the direction of the thickness of PCM layer. The effective heat capacity was parametrized as an asymmetrical Gaussian function: 302 ceff = c + c𝐿𝑓 exp {− (T − Tppc) 2 𝜎 } (2) where c is the specific heat capacity outside the phase change interval, c𝐿𝑓 is the hight coefficient defining the amount of the latent heat of fusion, Tppc stands for the peak phase change temperature and 𝜎 is the sharpness coefficient of ceff function defined as 𝜎 = { 𝜎𝑠, 𝑇 ≤ Tppc 𝜎𝑙 , 𝑇 > Tppc (3) where 𝜎𝑠 and 𝜎𝑙 are corresponding to sharpness’s of solid and liquid phase. The symmetrical case is equivalent to 𝜎 = 𝜎𝑠 = 𝜎𝑙 (see Figure 2, where 𝜎𝑠 < 𝜎𝑙). Figure 2: An asymmetrical Gaussian effective heat capacity curve Eq. 1 is solved by the energy balance approach following the forward time-explicit numerical scheme. As was previously mentioned, the numerical model consists of a series of 1-D sub-models with each model corresponding to a specific section of the CSM. The neighbouring sub-models were coupled by means of the model of the air flow creating quasi 2-D model of the unit. The numerical model of an air-PCM HX is based on authors’ previous publication (Charvát et al., 2014), where a similar approach was adapted and described in a much greater detail. The air-PCM HX of 100 CSMs was investigated both experimentally and numerically. The reader is also referred to (Stritih et al., 2018), where the thermal behaviour of a solar air heating system containing a solar collector and an air-PCM HX was studied. The division into a series of 1-D sub-models and the corresponding coupling with the air flow sub-model were closely explained in the study. Figure 3: The boundary conditions and the presimulated outlet air temperature 303 A paraffin-based Rubitherm RT series PCM was considered in the study. The corresponding material properties are summarised in Table 1 as the “presimulated” case. The boundary condition (the inlet air temperature, see Figure 3) was set to 55 °C and 15 °C during the melting and cooling phases. Both phases lasted 150 min each, which was enough time for the PCM to fully undergo the phase transition. The adiabatic boundary conditions were considered along the air-PCM HX outer surfaces. Consequently, all heat losses into the surroundings were neglected. The identified parameters were as follows: the shape defining parameters of ceff function c, c𝐿𝑓 , Tppc, 𝜎𝑠, 𝜎𝑙 (see Eq. 2), density 𝜌 and thermal conductivity 𝑘 of PCM, defining the vector of optimisation variables as p = (c, c𝐿𝑓 , Tppc, 𝜎𝑠, 𝜎𝑙, 𝜌, 𝑘). In order to assess the overall robustness and accuracy of the inverse problem optimisation, the presimulated scenario was constructed by the direct run of the numerical model for a given set of parameters 𝐩presim (shown in Table 1). The corresponding calculated outlet air temperature Tout,presim air (shown in Figure 3) was used instead of the experimentally obtained data. As a result, the existence of the global optima OF(p) = OF(𝐩presim) = 0 is guaranteed and the uncertainties tied to the measurement methods are avoided. As for the optimisation variables, in total 7 distinct parameter were considered. The objective function was defined as a sum of the squared differences between the presimulated Tout,presim air and calculated Tout,,sim air outlet air temperatures: OF(p) = ∑(Tout,,sim air − Tout,presim air ) 2 . tmax ∆t p=0 (4) 3. Results Several inverse problems were solved with the PSO method. The inverse problems mainly differed in the number of unknown variables (parameters) that were identified with PSO optimisation. A robust solution method to an inverse problem should be able to produce accurate results regardless of the number of unknown variables or the size of the state space. The inverse problems were constructed with the increasing number of unknown variables (parameters) as follows: 1. Case: Tppc – A single unknown parameter (the peak phase change temperature Tppc) was identified (considered as the optimisation variable). 2. Case: 𝜌, 𝑘 – Only parameters that are not directly tied to the shape of the effective heat capacity curve (𝜌 and 𝑘) were determined. 3. Case: ceff – The complementary case to the inverse problem 2, only ceff defining parameters (c, c𝐿𝑓 , Tppc, 𝜎𝑠, 𝜎𝑙) were used as the variables. 4. Complete Symmetrical – In total 6 unknown parameters were considered (c, c𝐿𝑓 , Tppc, 𝜎 , 𝜌, 𝑘), since the ceff is considered symmetrical (𝜎𝑠 = 𝜎𝑙 ). The presimulated scenario is asymmetric. Using the symmetrical ceff will be always an approximation as OF(p) ≠ 0. 5. Complete Asymmetrical – All parameters in 𝐩presim were unknown (identified). The effective heat capacity curve was asymmetrical (𝜎𝑠 ≠ 𝜎𝑙 ). 3.1 Optimisation The particle swarm optimisation (PSO) method was used as an optimisation solver to the inverse problems. The main advantage of PSO is mainly a good computational efficiency and a small amount of control parameters (Yang, 2014). The method is based on a random selection of individuals in the search space (particles) which update their position (in a feasible region) during each iteration according to their speed vector. The size of population Np was chosen as 𝑁𝑝 = 10 × 𝑁𝑣𝑎𝑟, where 𝑁𝑣𝑎𝑟 is the number of unknown variables in a given optimisation case, since (Piotrowski et al., 2020) have shown that 20 – 100 is usually the optimal population size for the lower dimension problems. As a termination criterion, the maximum number of 10 stale iterations (iterations with no improvement in the OF) was chosen, coupled with the objective value criterion OF < 0.05. The upper and lower optimisation bounds are show in Table 1. The accuracy of the optimisation solver was analysed with regard to the shape of the ceff function, the enthalpy of fusion 𝐿𝑓 and its percentage difference to the presimulated results and the value of the OF. The ceff- temperature curves are shown in Figure 4 with corresponding parameters being summarized in Table 1. There is no observable discrepancy in Case: Tppc and Case: No ceff (OF = 0.04 K2 in both cases) since the inverse problem had only 1 and 2 unknown variables. Possibly a method such as Nelder-Mead optimisation method can reach comparable accuracy with faster computational time (however the risk of being stuck in a local optimum 304 would increase). In the Case: ceff, Complete S. and Complete As., the objective was to find parameters of the ceff function. Since the presimulated results were obtained with the asymmetrical Gaussian ceff, there is an observable difference when the solution in the form of a symmetrical Gaussian function was sought (OF = 1.58 K2, see Figure 4 dashed cyan curve). The maximal difference in the enthalpy of fusion was 2 %. It was in case of the inverse problem where the effective heat capacity curve in the form of asymmetrical Gaussian curve was sought in the form of a symmetrical Gaussian curve. The effective heat capacity curves of real-life PCMs never have the shape of a symmetrical Gaussian curve; however, symmetrical Gaussian curves are often used as an approximation in the simulation models. The result of present study show that such approximation can provide quite accurate results. In order to assess the robustness of the optimisation method, the number of variables was increasing from 1 to 7 as shown in Table 1. The increasing number of unknown variables did not have a significant influence on the overall accuracy (see Fig. 4). However, the computational time increased with the main reason being the population size growth and the rise in the overall number of iterations Nit required to reach comparable accuracy. Table 1: The optimisation results, the optimisation bounds and the parameters of the PSO method Case Unit 𝑐 [J/(kg K)] c𝐿𝑓 [J/(kg K)] Tppc [°C] 𝜎𝑠 [K2] 𝜎𝑙 [K2] 𝜌 [kg/m3] 𝑘 [W/(mK)] 𝐿𝑓 [kJ/kg] OF [K2] Comp. time [s] Nit [-] Np [-] 𝐩presim 2,000 20,000 35 23 3 800 0.2 115.70 N/A N/A N/A N/A Case: Tppc N/A N/A 34.97 N/A N/A N/A N/A N/A 0.04 229 14 10 Case: 𝜌, 𝑘 N/A N/A N/A N/A N/A 800.1 0.19 N/A 0.04 491 16 20 Case: ceff 2,001.8 19,594 34.81 22 3.9 N/A N/A 115.74 0.03 1,697 27 40 Complete S. 1,999.2 19,172 33.35 12.05 12.05 791.31 0.19 117.96 1.58 9,632 235 60 Complete As. 2,002.1 20,443 35.48 26.51 1.46 802.06 0.21 115.17 0.09 9,117 227 70 Upper bound 500 5,000 10 0.5 0.5 200 0.01 N/A N/A N/A N/A N/A Lower bound 5,000 40,000 60 50 50 2,000 10 N/A N/A N/A N/A N/A Figure 4: The comparison of effective heat capacity curves 4. Conclusions The present study focused on the accuracy and robustness of the particle swarm optimisation method used for the identification of the material properties of PCM such as effective heat capacity, density and thermal conductivity. The outcomes of this study can be summarized as:  The identification procedure was developed by coupling of the numerical model of an air-PCM HX with the PSO method.  In order to analyse the procedure, 5 distinct inverse problems were defined varying in the number of unknown parameters (optimisation variables) from 1 to 7. 305  The overall accuracy was analysed using the OF value, shape defining parameters of ceff function and the percentage change in the enthalpy of fusion. The results have shown that the highest discrepancy was in case of the Complete S. scenario with OF = 1.58 K2 and 2 % discrepancy in the enthalpy of fusion. This was caused mainly by the fact, that the asymmetrical ceff function was used in the presimulated scenario. The symmetrical case was an approximation. All other cases did differ only by up to 0.5 % considering the enthalpy of fusion.  The computational time did increase with the higher number of optimisation variables due to the increasing size of the population and the number of iterations needed to reach comparable accuracy. Combining the shape defining parameters of the ceff function with density and thermal conductivity did lead to much higher computational times. As was shown in Case ceff, investigating ceff parameters separately did speed up the optimisation considerably.  The accuracy of the solutions to the inverse problems did not significantly decrease with the increasing number of unknown variables. Therefore, the robustness of the method has been demonstrated. Nomenclature 𝑐 – specific heat outside the phase transition, J/(kg K) c𝐿𝑓 – hight coefficient of the ceff, J/(kg K) ceff – effective heat capacity, J/(kg K) 𝑘 – thermal conductivity, W/(m K) 𝐿𝑓 – enthalpy of fusion, J/kg Nit – number of iterations Np – population size Nvar – number of variables OF – objective function, K2 𝐩 – vector of the optimisation parameters T – temperature of PCM, K Tppc – peak temperature of phase change, K Tout air – outlet air temperature, K x, y – cartesian coordinates 𝜌 - density of PCM, kg/m3 𝜎 – sharpness coefficient of the ceff, K2 Acknowledgements This work was supported by the project “Computer Simulations for Effective Low-Emission Energy Engineering”, No. CZ.02.1.01/0.0/0.0/16_026/0008392 and by the project “Sustainable Process Integration Laboratory – SPIL”, No. CZ.02.1.01/0.0/0.0/15_003/0000456, both funded by European Research Development Fund, Czech Republic Operational Programme Research, Development and Education, Priority 1: Strengthening capacity for quality research. References Cascone Y., Perino M., 2015. Estimation of the thermal properties of PCMS through inverse modelling. Energy Procedia, 78, 1714-1719. Charvát P., Klimeš L., Ostrý M., 2014. Numerical and experimental investigation of a pcm-based thermal storage unit for solar air systems. Energy and Buildings, 68, 488-497. Franquet E., Gibout S., Bédécarrats J.-P., Haillot D., Dumas J.-P., 2012. Inverse method for the identification of the enthalpy of phase change materials from calorimetry experiments. Thermochimica Acta, 546, 61-80. Gossard D., Karkri M., Al Maadeed M.A., Krupa I., 2015. A new experimental device and inverse method to characterize thermal properties of composite phase change materials. Composite Structures, 133, 1149- 1159. Grabo M., Weber D., Paul A., Klaus T., Bermpohl W., Krauter S., Kenig E.Y., 2019. Numerical Investigation of the Temperature Distribution in PCM-integrated Solar Modules. Chemical Engineering Transactions, 76, 895-900. Liu Z., Sun P., Xie M., Zhou Y., He Y., Zhang G., Chen D., Li S., Yan Z., Qin D., 2021. Multivariant optimization and sensitivity analysis of an experimental vertical earth-to-air heat exchanger system integrating phase change material with Taguchi method. Renewable Energy, 173, 401-414. Piotrowski A.P., Napiorkowski J.J., Piotrowska A.E., 2020. Population size in Particle Swarm Optimization. Swarm and Evolutionary Computation, 58, 100718. Stritih U., Charvát P., Koželj R., Klimeš L., Osterman E., Ostrý M., Butala V., 2018. PCM thermal energy storage in solar heating of ventilation air, experimental and numerical investigations. Sustainable Cities and Society, 37, 104-115. Yang X.-S., 2014. Chapter 2 - analysis of algorithms, Nature-Inspired Optimization Algorithms, Yang X.-S (Ed.), Elsevier, Oxford, UK, 23-44. Zálešák M., Klimeš L., Charvát P., 2020. Design Optimization of a Solar Air Collector Integrating a Phase Change Material. Chemical Engineering Transactions, 81, 211-216. 306