Mechanical reliability analysis of nanoencapsulated phase change materials combining Monte Carlo technique and the finite element method

Nanoencapsulated phase change materials (nePCMs) are one of the technologies currently under research for energy storage purposes. These nePCMs are composed of a phase change core surrounded by a shell which confines the core material when this one is in liquid phase. One of the problems experimentally encountered when applying thermal cycles to the nePCMs is that their shell fails mechanically and the thermal stresses arising may be one of the causes of this failure. In order to evaluate the impact of the uncertainties of material and geometrical parameters available for nePCMs, the present work presents a probabilistic numerical tool, which combines Monte Carlo techniques and a finite element thermomechanical model with phase change, to study two key magnitudes of nePCMs for energy storage applications of tin and aluminium nePCMs: the maximum Rankine’s equivalent stress and the energy density capability. Then, both uncertainty and sensitivity analyses are performed to determine the physical parameters that have the most significant influence on the maximum Rankine’s stress, which are found to be the melting temperature and the thermal expansion of the core. Finally, both a deterministic and a probabilistic failure criterion are considered to analyse its influence on the number of predicted failures, specially when dispersion on tensile strength measurements exists as well. Only 1.87% of tin nePCMs are expected to fail mechanically while aluminium ones are not likely to resist.


Introduction
The world demand of energy is estimated to increase by 26% and CO 2 associated emissions will continue rising by 10% by the year 2040 with respect to those registered in 2017 (International Energy Agen, 2019). In order to reduce the environmental problems, efforts are made by several scientific communities to foster the use of renewable energies, which exploit natural resources -unlimited on a human timescalewithout generating polluting emissions. The present work focuses exclusively on solar energy and more precisely, on its application to concentrated solar plants (CSP), since solar energy is the renewable energy presenting the major potential for exploitation of energy (International Energy Agen, 2011). Nevertheless, one of the main cons of solar energy is its intermittence since energy production depends largely on weather/climate conditions. Consequently, owing to the high demand of energy in the electricity market, gaps between generation and supply of energy cannot be tolerated to guarantee the correct operation of the electrical grid. One of the characteristics that makes CSP stand out among other renewable energy technologies is the possibility of incorporating Thermal Energy Storage (TES) systems to mitigate the previously mentioned generation gaps. TES systems appear to be a field of research on its own nowadays for energetic transition towards renewable energies within the context of green policies (Gil et al., 2010;Xu et al., 2015;Akhmetov et al., 2016;Mondragón et al., 2017).
The focus of the present work lays on one of the technologies under research for TES systems in concentrated solar plants: nanofluids (Choi and Eastman, 1995), which is the name received by the colloidal suspension of nanoparticles (less than 100 nm) in a base fluid to enhance its thermal properties. One of the advantages of these suspensions is that, due to the Brownian motion of the particles within the base fluid, nanofluids combine the good thermal properties of nanosolids with the transport properties of the fluid avoiding clogging or settling problems.
The thermal energy storage capability of the nanofluids is a key parameter of their performance for TES purposes. Although the enhancement of sensible heat storage capability of commonly used nanoparticles (metal oxides or carbon structures) dispersed in ionic liquids still generates controversy (Riazi et al., 2016), it has been recently demonstrated the possible enhancement of thermal energy storage of TES material with the addition of metallic nanoencapsulated Phase Change Materials (nePCMs) (Navarrete et al., , 2019Cingarapu et al., 2013Cingarapu et al., , 2015 due to the contribution of the latent heat. These nanoencapsulated particles are composed of a phase change core wrapped in a surrounding shell made of another material with a melting temperature higher than that of the core. The role of the shell is to confine the core material when melting occurs. One of the main features of metallic nePCMs for TES applications is that they contribute considerably to increase the overall energetic performance of the nanofluid thanks to the latent heat absorbed/released during the phase change of their core. Moreover, they can be self-encapsulated with an oxide shell formed by natural passivation during the sintering process that can interact with the base fluid increasing the sensible heat in ionic liquids. In this context, optimising the design of the nePCMs is of key importance given that they play a direct influence on the performance of the nanofluid. One of the main problems encountered when nePCMs undergo thermal cycles is that their shell may fail mechanically, as it has been reported in ). The thermal stresses arising in the shells of the nePCMs could be one of the reasons for this failure, as investigated in (Forner-Escrig et al., (2020)).
Despite the fact that synthesis of nePCMs is performed experimentally with the valuable knowledge of the experimental community, the current work presents a numerical tool to analyse the behaviour of nePCMs from a quantitative perspective to get a better understanding of the physical phenomena involved in the failure of the nanoparticle shells. Furthermore, experimental measurements of the different material properties of a substance are not exempt from uncertainties, which are intrinsic to the nature of the measurement process. For this reason, it is important to numerically check the influence and consequences of these uncertainties on the measurements of the material properties by using a probabilistic technique that combines both Monte Carlo (MC) (Kroese et al., 2011) and Finite Element (FE) methods. In the same vein, considerable dispersion exists normally in the mean size of the synthesised nePCMs  and despite possessing a spherical morphology, they are not perfect spheres. Therefore, this lack of symmetry could have an impact on the mechanical failure of the shell of the nePCMs, which needs to be quantitatively assessed. Consequently, the study of all the aforementioned uncertainties on the nePCMs highlights the need of considering probabilistic analyses to evaluate the performance of these nanoparticles for TES applications in a wide range of temperature such as cold TES, comfort applications or high-temperature thermal management systems.
Some examples of probabilistic and numerical analyses are reported in the literature to analyse the failure of polymeric nanoparticle composites (Hamdia et al., 2017), and to optimise the nanoparticle wet milling process (Hou et al., 2007). However, neither their scope of analysis is the same as the one intended in the present work (nePCMs) nor the number of input parameters of their probabilistic analysis is large enough to analyse the influence of different physical parameters on their aim of study. Notice that the present work considers a set of 22 physical parameters (18 material and 4 geometrical parameters) to account for its influence on the failure of the nePCMs and, consequently, the analytical computation of the probability of failure may be complex. Then, the need for numerical methods arises.
MC simulation is a class of algorithms that use statistical sampling techniques to obtain a probabilistic approximation to the solution of a model: the samples previously generated are the inputs of the model and its evaluations provide the probabilistic outputs or responses of the model (Aderibigbe, 2014). The evaluation of the model may be performed by using analytical solutions or numerical methods such as the FE (Pérez-Aparicio et al., 2012;Palma et al., 2013Palma et al., , 2018, finite difference (Li and Rankin, 2010), finite volume (Bhoi and Sarkar, 2020), etc.
On this ground, the aim of the current work is to study two key magnitudes of nePCMs for TES applications: the maximum equivalent stress and the energy density capability and, to achieve this goal, reliability, uncertainty and sensitivity analyses are developed. In particular, two different pairs of nePCM materials are considered: tin encapsulated in tin oxide and aluminium encapsulated in aluminium oxide (alumina). For this purpose, a probabilistic tool is obtained by combining both MC and FE methods. The former considers up to 22 random variables and the input samples are generated by using the Latin Hypercube Sampling (LHS) technique (McKay et al., 1979), while the latter is a thermomechanical with phase-change FE code previously developed by the authors of the present work (Forner-Escrig et al., 2020). From the uncertainty analysis, distributions of maximum Rankine's equivalent stress -failure criterion for the shell-and of energy density are obtained. Furthermore, the sensitivity analysis permits to assess which of the 22 random variables exerts more influence on both the failure of the nePCM and on its energy density performance. With regard to the reliability analysis, both a deterministic and a probabilistic failure criterion are considered to verify its influence on the number of predicted failures when uncertainty in tensile strength exists. In conclusion, the results of these analyses may contribute to identify what are the physical parameters to be considered in advance to synthesise mechanically reliable nePCMs.

Outline of the problem
This section briefly reviews the model of the nePCM retained for the probabilistic analyses henceforth.

Description of the nanoparticle
The numerical model considers a single three-dimensional ellipsoidal nePCM, as that sketched in Fig. 1 a), where the geometry is defined by the three semi-axes a, b and c of the outer ellipsoid and the shell thickness e shell . Notice that an ellipsoidal geometry for nePCMs is assumed in the present work for two reasons: i) to consider the influence of the geometrical uncertainty in the probabilistic analyses and ii) given that real nePCMs are not perfectly spherical. Fig. 1 b) depicts the mesh of an ellipsoidal nePCM and its appearance is quasi-spherical since the dispersion around the nominal values of the semi-axes is only of 5% (see Section 3.1 and Table 1).
The material properties listed in Table 1 are mass density ρ, specific heat capacity at constant pressure c p , thermal conductivity κ, Young's modulus E, Poisson's ratio ν, thermal expansion coefficient α, melting temperature T m and latent heat L. Subscripts s and l denote solid and liquid state, respectively. For the geometry, an order of magnitude of the size of spherical Sn and Al nanoparticles may be obtained from ) and (Sarathi et al., 2007;Sahu and Hiremath, 2017;Navarrete et al., 2019), respectively.
Concerning the boundary and initial conditions, the nePCM is mechanically fixed at its centre and subjected to an initial temperature T i . Then, an increasingly linear temperature is applied at the outer surface of the shell until a value of temperature T 0 (higher than the melting temperature of the core material) is reached. Concretely, T i = 323.15 (K), T 0 = 523.15 (K) for Sn@SnO 2 and T 0 = 973.15 (K) for Al@Al 2 O 3 nePCMs are considered.

Description of the finite element model
The evaluation of the nanoparticle model is performed by a threedimensional FE code developed by the authors. In particular, a thermodynamically consistent FE formulation was developed in (Forner-Escrig et al., (2020)) by considering thermomechanical coupling and three phase change approaches: enthalpy, heat source and heat capacity. This formulation was implemented in the research code FEAP (Taylor, 2013), which belongs to the University of California at Berkeley (USA).
Numerically, the FE formulation is monolithic and eight-noded with four degrees of freedom (dof) per node: where the first term in brackets represents the three components of the displacement field u = {u, v, w} and the last term is the temperature. These dofs are discretised by using standard shape functions of Lagrangian type. For the sake of brevity, only the tangent matrix is reported in the present work: ⎡ where K and M denote stiffness and mass matrices, respectively. Subscripts a and b refer to discretisation nodes while superscripts denote the dofs. As observed, only a one-way coupling is considered with the term K u T ab , which is due to the thermomechanical volumetric expansion. With regard to the time integration scheme, the scalar coefficient c 3 in (2) is associated to the Newmark− β algorithm. Finally, the enthalpy approach is used in the current work.

Probabilisitic analysis
The present work uses a probabilistic numerical tool that combines the thermomechanical FE method with the MC technique. As shown in

Distribution functions and sample generation
Firstly, the j random variables ξ j -commonly called parameters-of the model must be identified and quantified in accordance with experimental observations. In the present work, all the material and geometrical properties (j = 1, …, 22) of the nePCM are considered as random variables; their nominal values and standard deviation (uncertainty) of the random variables are shown in Table 1.
Secondly, neither the distribution function of the measurements of these variables nor most of uncertainties in the measurement process are available in the literature. However, the experimental values of physical properties are normally obtained as the average of the different measurements and, therefore, according to the central limit theorem, these values can be assumed to be normally distributed even when the original samples themselves do not obey a normal distribution.  According to the orders of magnitude of dispersion in measurements reported in (Perry et al., (2008)), a standard deviation σ of 5% around the nominal value (mean) of each material parameter is considered as a first and good estimation. With regard to geometrical parameters and in agreement with the experimental characterisations available in (Navarrete et al., , 2019, the size of nePCMs is assumed to be log-normally distributed. In particular, the mean size and standard deviation of the log-normal distribution for the generation of the nominal value of the outer ellipsoid size are 40 [nm] and 40% around the mean value, respectively. Since experimentally it is observed that nanoparticles are not perfect spheres but quasi-spherical, an ellipsoidal geometry is defined with a standard deviation of 5% around the previously generated log-normal values of the nanoparticle size. Notice that the shell thickness is uniform around the core of a single nePCM but its nominal value can vary for the different nePCMs within a 5% of uncertainty by following a normal distribution. Finally, the LHS technique (McKay et al., 1979) is used to generate the random samples. LHS is a statistical technique to generate random numbers which consists of dividing the cumulative density function of each random variable into N equal partitions and choosing a random data point in each partition of each random variable. Then, the samples of each random variable are combined in order to create the set of random input vectors of the model. This technique is used in the present work for two reasons: i) it reduces the computational time (Olsson et al., 2003) with respect to the time required for standard random sampling, and ii) the random sample generated is more representative of the variability of the random variables than standard random generations.

Uncertainty and sensitivity analyses
Consider a generic model M represented by the mathematical expression: where ϕ i represents the i outputs and ξ j the j inputs. On the one hand, the main goal of the UA is to determine the uncertainty in the model output when the uncertainties in the input are known. Concretely, UA allows to calculate the probability distribution of the output and their scalar parameters: mean and standard deviation.
On the other hand, the objective of the SA is to determine the relationship among the uncertainties in input and output variables, namely, SA permits to identify which input variables exert more influence on the response of the model. Although there are different techniques to develop a SA (Saltelli et al., 2008), multiple linear regression is adopted in the present work. According to (Montgomery and Runger, (2003); Palma et al., 2009), in multiple linear regression, each output ϕ i may be approximated by: where N ξ represents the number of random variables (or inputs) of the model, θ j are the regression coefficients that relate each input ξ j with each output ϕ i , θ 0 is the intercept term of the linear regression and Ψ represents a random error term. Regression coefficients θ 0 and θ j are determined by least-square computation from (4) and, despite the fact that they provide information on the relationship between the inputs and outputs of the model, according to (Saltelli et al., (2008)), regression coefficients are not an appropriate sensitivity measure. Instead, another sensitivity measure known as Standardised Regression Coefficients (SRC) is preferred and they are defined as: where σ ξ j and σ ϕ i denote the standard deviations of input and output, respectively. SRC are recommended measures for SA since they account not only for raw regression coefficients but also for standard deviations of both inputs and outputs and then provide a normalised measure of the importance of the input parameters on the response of the model (Saltelli et al., 2008). Thus, SRC are the sensitivity measures considered through the present work.

Probability of failure
The POF, which is a statistical indicator of the frequency of occurrence of a considered failure event, is mathematically defined as (Melchers and Beck, 2018): where ξ j , G(ξ j ) and f ξ j (Ξ) represent the vector of input random variables, a limit state function and the joint probability density function of the input random variables, respectively. Therefore, G(ξ j ) ≤ 0 refers to the region where the limit state violation occurs, i.e. the failure region. In general, equation (6) cannot be evaluated analytically, except for some special cases, but the POF can still be determined numerically. One of the existing techniques to compute this POF numerically is MC techniques and according to (Melchers and Beck, (2018)), for these techniques, the POF can be evaluated as: where n[G(Ξ j ≤ 0)] is the number of trials n for which the limit state function is not satisfied, i.e. G(Ξ j ≤ 0). In (7), Ξ j is used to represent a sample value of each random variable.
As well known, the number of total trials N determines the desired accuracy of the POF, which may be measured by the Coefficient of Variation (CV) of the POF (Smarslok et al., 2010): Then, the larger N, the most accurate POF is but, in contrast, choosing large values of N increases considerably the computational cost and, therefore, an agreement between accuracy and computation time must be reached. The value of N used in the present work and the time needed to run these simulations is detailed in Section 4.1.

Results
This section summarises the results obtained from the probabilistic tool for both Sn@SnO 2 and Al@Al 2 O 3 nePCMs.

Uncertainty analysis
Two UA for two nePCMs -Sn@SnO 2 and Al@Al 2 O 3 -are conducted in order to obtain two outputs of the model (3): i) The maximum Rankine's equivalent stress σ R distribution on the nePCM shell: ϕ 1 = σ R , see Fig. 3.
ii) The energy density E d of the nePCMs: ϕ 2 = E d , see Fig. 4.
Notice that these outputs are selected among others since they are key magnitudes in the design and optimisation of nePCMs for TES applications.
With regard to the equivalent stress distribution and in agreement with (Forner-Escrig et al., 2020), Rankine's criterion is used to predict the mechanical failure of the nePCM, which usually occurs at the shell. Fig. 3 shows the histogram and the statistical scalar values (mean μ and standard deviation σ) after running the UA. results, their tendency appear to be logical since larger values of stress are more likely to occur in Al@Al 2 O 3 . Notice that the melting temperature of the Al core is considerably higher than that of Sn; consequently, larger thermal stresses are expected to arise in the Al@Al 2 O 3 nePCMs to reach the liquid state of the core.
In order to measure the accuracy of the predicted σ R , a bilateral Confidence Interval (CI) around the mean μ is defined as: where α represents the risk -also called significance level-of the CI and t N− 1 is the Student's t-distribution with N − 1 degrees of freedom. The Standard Error (SE) of μ is defined as: The CI for the equivalent stress at a significance level of α = 5% is Regarding energy density E d , which is a relevant magnitude for energy storage purposes, the metric used in the present work for its quantification is defined as: where V core and V total denote the core and total volume of the nePCM, respectively. Fig. 4 represents the histogram and scalar values of E d for Sn@SnO 2 Fig. 3. Histogram of Rankine's equivalent stress for Sn@SnO 2 nePCMs (left) and Al@Al 2 O 3 nePCMs (right) with their respective mean value μ and standard deviation σ. concluded that the energetic performance of Al nePCMs is better than that of Sn ones. Two numerical verifications are performed: statistical and convergence tests. On the one hand, the Jarque-Bera test (Jarque and Bera, 1987) is used to determine if the null hypothesis -a data sample comes from a normal distribution at a certain significance level-is accepted.
The Jarque-Bera test is applied to both σ R and E d and it is found that both distributions and both pairs of core@shell materials are normally distributed at a significance level of 1%. Notice that the linearity assumption is not obvious since the present probabilistic analysis combines normal and log-normal distributions.
On the other hand, the sample size N to guarantee a proper convergence of the MC is performed by a trial and error numerical test since obtaining a general analytic expression to pre-calculate N for convergence of MC simulations is difficult or even impossible for complex models. The numerical test consists of representing the average evolution of μ and σ versus the number of iterations, as observed in   is approximately 14 h for a PC with a processor Intel® Core™i7-950 and 16 GB of RAM.

Sensitivity analysis
A multiple linear regression is performed to develop the SA in order to calculate the influence of the random variables ξ j on the σ R for both Sn@SnO 2 and Al@Al 2 O 3 nePCMs. This influence is shown in Fig. 6 by a bar diagram for the different absolute values of the SRC, see notation in Table 1.
Firstly, the SRC values obtained for each nePCM are not the same given that the material properties are different. However, the same tendency may be observed in both of bar diagrams.
Secondly, the most relevant random values appear to be Θ 8 , Θ 9 and Θ 10 , which represent the Poisson's ratio of the core, the thermal expansion coefficient of the core and its melting temperature, respectively. These results seem to be in good agreement with physical intuition: the latter plays a direct role in the maximum value of stress reached until the melting of the core, and the thermal expansion coefficient and the Poisson's ratio are parameters that mechanically govern the volumetric changes in nePCMs.
Thirdly, the next set of parameters having more influence on the failure are the geometrical ones: Θ 19 , Θ 20 , Θ 21 and Θ 22 . The three first parameters stand for the semi-axes of the ellipsoidal nePCM and the last one represents the shell thickness. The difference in values between Sn@SnO 2 and Al@Al 2 O 3 nePCMs for the SRC of the three semi-axes in Fig. 6 is not very significant and could be influenced by randomness in the generation of the random samples. Finally and since multiple linear regression is applied for this sensitivity analysis, the validity of the hypothesis of linearity must hold true for the results to be acceptable. In order to verify the validity of the SA, the coefficients of determination of the linear regression are R 2 Sn = 0.93 and R 2 Al = 0.92 for the Sn@SnO 2 and Al@Al 2 O 3 , respectively. According to (Pan et al., (2011);Saltelli et al., (2006)), as long as the coefficient of determination verifies R 2 ≥ 0.7, the linearity hypothesis is satisfied and thus, the present SA is valid as well.
Notice that the random variables ξ j exerting a major influence on E d can be directly determined from its definition in Equation (11) without need of performing a SA. Therefore, the latent heat of the core, the mass density of the liquid core and the geometrical parameters are the variables that have a more relevant impact on E d .  Table 1.

Probability of failure of nePCMs
In this section, the POF of Sn@SnO 2 and Al@Al 2 O 3 nePCMs is calculated by using both deterministic and probabilistic failure criteria: • The deterministic criterion consists of using a tensile strength σ t value from literature and studying its position with respect to the σ R distribution obtained from the probabilistic simulations. • The probabilistic criterion is grounded on the uncertainty in the measure of σ t , resulting in a tensile strength distribution. According to (8), by considering a probabilistic criterion to compute the POF, the accuracy of the MC simulation is increased by a factor ̅̅̅ ̅ N √ -a factor of 10 in the present simulations. Notice that this motivates the need for considering a probabilistic failure criterion.
Concerning Sn@SnO 2 nePCMs, σ t = 803 [MPa] is assumed for SnO 2 (Nam et al., 2017) for the deterministic case. According to (Nam et al., (2017)), the tensile strength of SnO 2 may suffer from fluctuations due to its porosity, which is a hard to control parameter when synthesising nePCMs. Consequently, σ t is assumed to be normally distributed with a 20% of dispersion around σ t = 803 [MPa] for the probabilistic criterion.
This percentage of 20% is considered as a first approximation from the measurement dispersion existing in the characterisation of tensile strength of metals at the macroscale and by assuming that a considerable error may be induced by the different behaviour of this mechanical property between the macro and the nanoscale. Fig. 7 shows the σ R distribution for both deterministic (top) and probabilistic (bottom) criteria. For the deterministic criterion, σ R falls below the deterministic σ t and, consequently, the POF of is 0%: any nePCM is expected to fail.
On the contrary and for the probabilistic criterion, it can be observed that there is an overlap between the Gaussian curve of the σ R distribution and that of the probabilistic σ t . This overlapping region represents the area in the graph where the nePCMs are likely to fail mechanically. In this case, the POF becomes 1.87% and this means that it is advisable to account for these type of criterion, especially when dispersion in tensile strength exists. It is important to highlight that the present results are in accordance with the experimental observations reported in ), given that only small samples of Sn@SnO 2 nePCMs were verified to fail mechanically due to thermal stresses.
With regard to Al@Al 2 O 3 nePCMs, to the best of the author's knowledge, experimental data on the mechanical failure of Al@Al 2 O 3 nePCMs is not available in the literature. In the present work, σ t = 275.9 [MPa] (Shackelford et al., 2015) with a Gaussian dispersion of 20% is considered, as in the previous case. Fig. 8 shows the distribution of the σ R with the deterministic (top) and the probabilistic (bottom) failure criteria. For this material, the POF predicted for both deterministic and probabilistic criteria is 100% since the distribution of σ R obtained is far above the tensile strength.
In conclusion, despite the fact that Al@Al 2 O 3 nePCMs possess larger energy density storage capabilities than Sn@SnO 2 nePCMs, the first ones are likely to fail mechanically when they undergo heating cycles.

Conclusions
This article has presented a numerical tool, which combines Monte Carlo techniques and a thermomechanical finite element with phase change, in order to conduct probabilistic analyses -uncertainty, sensitivity and reliability-in nanoencapsulated phase change materials (nePCMs). In particular, experimental uncertainties are taken into account to obtain the mechanical probability of failure and this mechanical failure is one of the problems experimentally encountered when the nePCMs undergo thermal processes. In addition, the sensitivity analysis has allowed to quantify the parameters that exert the most significant influence on the mechanical failure of nePCMs. Consequently, these relevant parameters should be closely controlled in the synthesis of nanoparticles.
Specifically, the present work has considered 22 random parameters − 18 material properties and 4 geometrical parameters-and has studied the evolution of two variables -Rankine's equivalent stress and energy density-for two nePCMs: Sn@SnO 2 and Al@Al 2 O 3 . Firstly and from the sensitivity analysis, it has been concluded that the melting temperature, thermal expansion coefficient and Poisson's ratio exert the more relevant influence on the mechanical failure of the nanoparticle shell. In turn, the most influential parameters on the energy density capability of nePCMs are the latent heat of the core, the mass density of the liquid core and the geometrical parameters of the nePCM. Secondly, it has been concluded that the probabilistic failure is 1.87% and 100% for Sn@SnO 2 and Al@Al 2 O 3 nePCMs, respectively. The results of the former are verified to be in good agreement with experimental observations while, to the author's knowledge, no data is available for the latter ones.
To sum up, the numerical probabilistic tool allows to estimate the probability of failure of the nePCMs and to determine the parameters having a bigger influence on the failure and energy storage of nePCMs for their use in TES systems. The scope of the present numerical tool can be expanded to predict the mechanical failure of any micro and nanoencapsulated PCMs made of different pairs of core@shell materials as long as their material properties and measurement dispersion are known.
Finally, the present numerical probabilistic tool could be used to complement experimental research and to reduce the number of experiments to be conducted to optimise the selection of a pair of core and shell nePCMs by maximising the energy density and by minimising the probability of the failure of the shell.