Basic Analysis of Uncertainty Sources in the CFD Simulation of a Shell-and-Tube Latent Thermal Energy Storage Unit

: Computational Fluid Dynamics (CFD) simulations are increasingly employed in the development of latent thermal energy storage units. Yet there are often strong deviations between the experiments and numerical results. To unveil the sources of the deviations for the CFD simulation of a vertical shell-and-tube latent thermal energy storage unit, a basic analysis of di ﬀ erent uncertainties is undertaken in this paper. Consequently, the e ﬀ ect of a variation of 10 material properties, six initial and boundary conditions, as well as a displacement of the temperature measuring points in the simulation, are examined. The results depict that the inﬂuence of the substance data depend on the output variable under consideration. Beside material properties, which have almost no inﬂuence, there are some properties that inﬂuence the power and the global liquid phase fraction over time, and a third group, which also has an inﬂuence on the mean power. Partly in contrast to results found in literature, the highest inﬂuence on the mean power occurs for the heat losses (which are varied in an on / o ﬀ manner), the density, and the melting enthalpy (both varied by ± 10%).


Introduction
The increasing share of renewable energies raises the temporal mismatch between energy production and consumption. As a consequence, energy storage systems attract more and more attention due to their ability to store and deliver energy when it is needed. Thermal energy storage systems are a cost effective option, not only for storing thermal energy, but they are also promising for demand side management, as well as for pumped thermal electricity storage [1]. If thermal energy needs to be stored in a narrow temperature range, or the energy storage density is of major importance, latent thermal energy storage systems (LTESS) are the right choice [2]. The high storage density of these storage units is a result of exploiting the latent heat of a phase change. Usually, solid/liquid phase changes of the so-called phase change materials (PCM) are chosen for this purpose. Analyzing, designing, and optimizing LTESS is more and more supported by numerical simulations. Thus far, however, major uncertainties have existed in these simulations, mainly due to a large variety of material properties reported in the literature [3], guessed constants in the numerical models, and a lack of highly precise validation experiments for most applications. In addition, it seems that these uncertainties are interpolating the pressure. There were no differences found in the results achieved with the SIMPLE or the PISO scheme, but the PRESTO! approach showed a better agreement with the experiments than the body-force-weighted approach. Lastly, it should be noted that independency studies for the mesh and the time step are commonly reported within papers on LTESS simulation, and a few authors compare two-dimensional (2D) with three-dimensional (3D) simulations showing that distinctive differences in the flow patterns may appear between both approaches [22,23].
Analyses regarding the material properties mostly only focused on one or a few specific aspects. From analyses with simple analytical models, Günter et al. [24] drew some basic conclusions regarding the effects of uncertain material data. The effect of differential scanning calorimetry (DSC) results that were achieved with different heating rates were investigated numerically [25,26], and the influence of the width of the melting range on the numerical results was studied by several authors [9,10,15,21,27,28]. Interestingly, it is reported that a larger melting range may either lead to a faster phase change [9,21] or a slower one [9]. Galione et al. [22] compared simulation results achieved with temperature dependent material properties to results that were obtained with constant material properties. It seemed that next to the density change between solid and liquid octadecane, the variation of the thermal conductivity of the liquid PCM has the largest influence on the melting rate. Soni et al. [29] checked the impact of assuming identical material properties in the solid and the liquid state compared to differing values for some of the material properties. The authors claim that the model with different material properties in the solid and liquid phase has better potential to reduce the discrepancy between simulation and experiment than the model with constant material properties. The effect of the uncertainty of the material properties data was analyzed in detail with purely diffusive phase change models by Dolado et al. [27], Zsembinszki et al. [30], and Mazo et al. [28]. In all three studies, the melting temperature is the most important, or one of the most important, material properties. Interestingly, the detailed uncertainty analyses performed in Dolado et al. [27] and Mazo et al. [28] allow to propose a certain set of allowed uncertainties in the material properties to achieve a predefined maximum error in the numerical results. Possible variances in the start and boundary conditions, as well as in the geometric data, are investigated by Dolado et al. [27] and Zsembinszki et al. [30]. Combined with the detailed analyses for the material properties, this allows for a fair comparison of the numerical results with experimental results as the overlap of both results can be taken as a measure for the agreement.
In summary, there are still large uncertainties present when it comes to the simulation of an LTESS, and even a basic validation of models for melting processes is still subject to distinct uncertainties. On the other hand, it is common practice to fit model constants, such as the Darcy constant to experimental results, and/or to choose material properties that fit the experimental results best. However, errors due to uncertain material properties, too coarse meshes and incorrect model assumptions as well as uncertainties in the experiment may have been hidden by a new error, which gives the illusion of a highly precise simulation. As stated above, there are already several studies on the influence of different enthalpy methods or switch-off techniques. However, to the authors' knowledge, there is no Computational Fluid Dynamics (CFD) study on LTESS so far that takes into account the uncertainties of the material properties, of the initial and boundary conditions, as well as of the experiment. Therefore, the aim of this work is to investigate these uncertainties with a basic sensitivity analysis performed with a CFD model describing an experiment from literature [31].

Set-Up
A vertical shell-and-tube type latent thermal energy storage unit is used as the set-up. The storage unit consists of two pipes. The outer pipe is made of Plexiglas and has an inner and an outer diameter of 44 mm and 50 mm, respectively. The inner and outer diameter of the inner pipe (made of stainless steel) are 10 mm and 15 mm, respectively. The height of the storage unit is 400 mm and a detailed description of the storage system is given in Longeon et al. [31]. A basic scheme of the storage system     For the heat transfer fluid (HTF) water and the steel pipe, the material properties available in ANSYS Fluent are used. The heat conductivity of the Plexiglas is set to λ Plexi = 0.19 W mK and the material properties of the PCM RT 35-nowadays called RT 35 HC-(Rubitherm Technologies GmbH, Berlin, Germany) can be found in Table 2.
The viscosity was measured with the IMETER Mess System and approximated by a linear fit.
For c s a high value appears because already below T s increased values for the heat capacity occur. The value of 5 kJ kg·K was measured at 30 • C by Rösler [32].

Boundary and Initial Conditions
At the beginning of every charging and discharging process, the storage system has a uniform temperature. The storage system is charged from the top and discharged from the bottom. After charging and discharging, the phase change process has already been finished in the reference simulation. The initial and boundary conditions of the reference simulation can be seen in Table 3. According to the publication of Longenon et al. [31], it is assumed in the reference case that the outer wall is adiabatic. Throughout the variations, t char as well as t dis are kept constant, and T start , T in and u in are varied by a fixed value. Additionally, T in is integrated once as a time curve, for u in a parabolic velocity profile is implemented once at the inlet and the heat losses are considered once. The time curves for T in are fitted to the measured inlet temperatures leading to for charging and for discharging. This results in maximum inlet temperature variations of about 1.7 K and 1.1 K for charging and discharging, respectively. To estimate the influence of the heat losses, a heat loss coefficient is defined, which includes the thermal resistance of the Plexiglas wall and the natural convection at the outside.
The natural convection at the outer surface A is calculated by means of laminar natural convection for a vertical cylinder [33]. Assuming T W,i = (T s + T l )/2 = 35.25°C on average and T amb = 18.5°C one gets for the overall heat flux coefficient U loss = 3.87 W m 2 K .

Numerical Model
The numerical model was set up in ANSYS Fluent 15.0 (ANSYS, Inc., Canonsburg, PA, US). The simulation domain consists of three different areas-the HTF in the inner pipe, the wall of the inner pipe, and the PCM in the annular gap. Within the HTF and the PCM, a laminar flow is assumed and the density is set constant and the Boussinesq approximation is used for the PCM. An enthalpy porosity method is applied [8] for the PCM, which implies two source terms. One in the energy equation to account for the latent enthalpy and another one, a Darcy-type source term, in the momentum equation to switch off the velocity in the solid region. The continuity equation with constant density is: and the momentum equation of the PCM becomes where the Boussinesq approximation is implemented by and the Darcy term is written as The energy equation of the PCM is In Equations (8) and (9) α F describes the liquid fraction, which is defined as and limited to values from 0 to 1. The update technique is described in detail in literature [34].

Mesh and Solver
The simulation domain was simplified. Only two dimensions are considered and an axis symmetry was used (see Figure 2). The mesh is orthogonal, has a maximum growth rate of 1.2 in y-direction in the HTF domain and cells with constant size in the domains of the pipe wall and the PCM. Since the mesh consists of 100 cells in x-direction and 122 cells in y-direction, it is ensured that a much finer resolution is generated in the expected direction of the largest gradients of velocity, temperature, enthalpy, etc. All relevant solver settings of the final simulations can be found in Table A1 in the Appendix A.
The mesh, the pressure-velocity coupling, the discretization of the convective terms in the energy, and momentum equation, as well as the time step will be analyzed in Section 3.1. Higher order methods for time discretization showed an unstable behavior. The tolerance for the residuals, the number of maximum iterations and the pressure interpolation scheme will not be varied. The mesh, the pressure-velocity coupling, the discretization of the convective terms in the energy, and momentum equation, as well as the time step will be analyzed in Section 3.1. Higher order methods for time discretization showed an unstable behavior. The tolerance for the residuals, the number of maximum iterations and the pressure interpolation scheme will not be varied.

Sensitivity Analysis
Within this work, a basic sensitivity analysis is performed where all variations are made with respect to one reference simulation. The parameters varied within the sensitivity analysis are subdivided into material properties as well as boundary and initial conditions (see Table 4). It should be mentioned that all varied input parameters are considered as uncertainties of the simulation. By doing so all of the uncertainties-except for the measurement position-can be compared with each other. In addition, the thermocouple positions within the PCM were varied within the simulations by +1 mm/−1 mm in their radial and vertical position. The variation values were chosen uniformly, as this allows easier transferability to other studies. Although here the variations refer to uncertainties, they can also give a first indication of the material properties for which a neglect of the temperature dependence could have a large impact.

Sensitivity Analysis
Within this work, a basic sensitivity analysis is performed where all variations are made with respect to one reference simulation. The parameters varied within the sensitivity analysis are subdivided into material properties as well as boundary and initial conditions (see Table 4). It should be mentioned that all varied input parameters are considered as uncertainties of the simulation. By doing so all of the uncertainties-except for the measurement position-can be compared with each other. In addition, the thermocouple positions within the PCM were varied within the simulations by +1 mm/−1 mm in their radial and vertical position. The variation values were chosen uniformly, as this allows easier transferability to other studies. Although here the variations refer to uncertainties, they can also give a first indication of the material properties for which a neglect of the temperature dependence could have a large impact.
For the evaluation, the simulation results of the variations are compared to the reference results without variation. The values that are compared are the absolute mean and the maximum of the absolute deviation of the power and the maximum of the absolute deviation of the global liquid fraction. The maximum deviations are calculated by and the mean deviations are calculated by Appl. Sci. 2020, 10, 6723 The time curves of the compared variables are Y re f and Y var for the reference case and the variation, respectively. The variable ∆Y is either the maximum of the global liquid fraction of the reference case (i.e., α F = 1) or the mean and maximum of the power of the reference case. It should be noted that with a constant time a change in average power (Equation (12)) corresponds to a change in the total heat transferred. Table 4. Varied material properties and initial and boundary conditions.

Validation and Independency Studies
All independency studies and the validation were performed for the charging process. The influence of the mesh was studied by three times refining an initial mesh of 12,200 cells with a factor of 1.2 in both spatial directions and by three times coarsening with a factor of 1.5 in both spatial directions. The maximum difference in the global liquid fraction (at 2000, 4000, and 6000 s) between the 12,200 cells mesh to the next finer one was only 0.63%, but the simulation time increased by about 40 h. Therefore, for all other simulations, a mesh with 12,200 cells was taken.
For solving the pressure-velocity coupling, the SIMPLE and the PISO scheme as well as a PISO scheme with extrapolation-here called PISO-extra-were investigated. No difference was found in the liquid fraction for the different schemes, but the PISO-extra scheme was the fastest one in terms of computational time and almost 46% faster than the slowest one (SIMPLE). All following simulations were performed with the PISO-extra scheme.
The influence of the discretization of the convective fluxes within the momentum and energy equation was studied as well. Three methods were analyzed (in an ascending sequence according to their order): Power-Law, Second-Order-Upwind, and Quadratic Upstream Interpolation for Convective Kinematics (QUICK). The maximum difference in the global liquid fraction at 2000, 4000, and 6000 s was about 4.5% between the Power-Law and the Second-Order-Upwind scheme, and about 0.5% between the Second-Order-Upwind and the QUICK scheme. Accordingly, the QUICK scheme was used for all subsequent simulations.
The time step was first calculated according to the ANSYS Fluent Guide [35] for natural convection to 0.0213 s. Subsequently, the time step was gradually increased and the results were compared. The simulation ran stable as long as a time step ≤ 0.2 s was chosen and there were no visible differences in the results for all stable simulations. In addition, a reduction of the simulation time of about 85% occurred by increasing the time step from 0.0213 to 0.2 s. For this reason, a time step of 0.2 s is used for all further simulations.

of 24
Finally C I was varied from 10 5 to 10 8 . It was seen that the results are not significantly different for a C I of 10 7 and 10 8 . Therefore, a value of 10 7 was chosen.
To validate the model, numerical results were compared to melting experiments described in Longeon et al. [31]. As depicted in Figure 3, the height of the solid PCM in contact with the outer tube is in good agreement between experiment and simulation. Due to the opacity of the solid PCM, the shape of the melting front cannot be analyzed optically. A comparison of the temperature curves over time, though, shows some significant deviations. This is illustrated in Figure 4 for temperature curves from measuring points with an identical height, and in Figure 5 for temperature curves from measuring points with an identical radius. Besides, the absolute mean discharging power up to 7000 s is 10.6% higher in the simulation compared to the experiment and the charging power is 29.9% smaller in the simulation. In order to better assess the reason for these deviations, a sensitivity analysis is performed in the following chapter.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 9 of 25 Finally was varied from 10 to 10 . It was seen that the results are not significantly different for a of 10 and 10 . Therefore, a value of 10 was chosen. To validate the model, numerical results were compared to melting experiments described in Longeon et al. [31]. As depicted in Figure 3, the height of the solid PCM in contact with the outer tube is in good agreement between experiment and simulation. Due to the opacity of the solid PCM, the shape of the melting front cannot be analyzed optically. A comparison of the temperature curves over time, though, shows some significant deviations. This is illustrated in Figure 4 for temperature curves from measuring points with an identical height, and in Figure 5 for temperature curves from measuring points with an identical radius. Besides, the absolute mean discharging power up to 7000 s is 10.6% higher in the simulation compared to the experiment and the charging power is 29.9% smaller in the simulation. In order to better assess the reason for these deviations, a sensitivity analysis is performed in the following chapter.   Finally was varied from 10 to 10 . It was seen that the results are not significantly different for a of 10 and 10 . Therefore, a value of 10 was chosen. To validate the model, numerical results were compared to melting experiments described in Longeon et al. [31]. As depicted in Figure 3, the height of the solid PCM in contact with the outer tube is in good agreement between experiment and simulation. Due to the opacity of the solid PCM, the shape of the melting front cannot be analyzed optically. A comparison of the temperature curves over time, though, shows some significant deviations. This is illustrated in Figure 4 for temperature curves from measuring points with an identical height, and in Figure 5 for temperature curves from measuring points with an identical radius. Besides, the absolute mean discharging power up to 7000 s is 10.6% higher in the simulation compared to the experiment and the charging power is 29.9% smaller in the simulation. In order to better assess the reason for these deviations, a sensitivity analysis is performed in the following chapter.

Temperature Measuring Position
The influence of the temperature measuring position is studied first. A variation of the measurement position in the numerical model by ±1 mm in radial position results in a difference in the PCM temperature of up to 5.2 K and a difference in the temporal mean of the temperature of up to 1.0 K. Varying the position in vertical direction by ±1 mm leads to much lower differences of a maximum deviation of less than 0.5 K and only about 0.03 K temporal mean deviation. Table A2 in the Appendix A shows the maximum and mean values as well as the standard deviation of the temperature differences introduced due to the measurement point variation. The outcome that a variation of the measurement position in radial direction has a much higher effect than a vertical variation can be explained by the almost vertical melting front (see Figure 3). The isotherms within the solid and the liquid PCM follow, to a large extent, the shape of the melt front and, therefore, a small and a large temperature gradient result in vertical and radial directions, respectively.

Material Properties
The outcome of the variation of the material properties is presented hereafter and an overview of the results including the sign of the resulting changes in the mean power can be found in the Appendix A in Table A3. Figure 6 shows the absolute relative variation of the mean power for charging and discharging with varied material properties. The variation of the density has the highest influence within the material properties of about 7.1% to 9.2%, followed by the melting enthalpy with 5.2% to 6.4%. The influence of the melting temperature, heat capacities and heat conductivity is always below 2.0%, but still noticeable. The influence is higher during charging compared to discharging for density, melting enthalpy as well as heat capacities and lower for heat conductivity and melting temperature. A variation of , , the Darcy constant and the mushy zone width − by ±10% and ±1 K, respectively, has almost no influence on the mean power for the chosen conditions. In general, the impact is somewhat higher when lowering the parameter than when increasing it.

Temperature Measuring Position
The influence of the temperature measuring position is studied first. A variation of the measurement position in the numerical model by ±1 mm in radial position results in a difference in the PCM temperature of up to 5.2 K and a difference in the temporal mean of the temperature of up to 1.0 K. Varying the position in vertical direction by ±1 mm leads to much lower differences of a maximum deviation of less than 0.5 K and only about 0.03 K temporal mean deviation. Table A2 in the Appendix A shows the maximum and mean values as well as the standard deviation of the temperature differences introduced due to the measurement point variation. The outcome that a variation of the measurement position in radial direction has a much higher effect than a vertical variation can be explained by the almost vertical melting front (see Figure 3). The isotherms within the solid and the liquid PCM follow, to a large extent, the shape of the melt front and, therefore, a small and a large temperature gradient result in vertical and radial directions, respectively.

Material Properties
The outcome of the variation of the material properties is presented hereafter and an overview of the results including the sign of the resulting changes in the mean power can be found in the Appendix A in Table A3. Figure 6 shows the absolute relative variation of the mean power for charging and discharging with varied material properties. The variation of the density has the highest influence within the material properties of about 7.1% to 9.2%, followed by the melting enthalpy with 5.2% to 6.4%. The influence of the melting temperature, heat capacities and heat conductivity is always below 2.0%, but still noticeable. The influence is higher during charging compared to discharging for density, melting enthalpy as well as heat capacities and lower for heat conductivity and melting temperature. A variation of η, β, the Darcy constant C I and the mushy zone width T l − T s by ±10% and ±1 K, respectively, has almost no influence on the mean power for the chosen conditions. In general, the impact is somewhat higher when lowering the parameter than when increasing it. discharging power over time is depicted in Figure 7. The highest influencing parameter is density with 13.3% to 18.4%, followed by melting enthalpy with 10.8% to 13.6%, heat conductivity with 7.9% to 9.8% and melting temperature with 5.5% to 7.8%. All other material properties have always an influence of less than 6.0% and has almost no impact. For , , , and especially c the influence is higher during discharging than during charging and for it is vice versa. Again, the impact is, in general, somewhat higher when lowering the parameters than when increasing them.   The maximal absolute relative influence of a variation in material properties on the charging and discharging power over time is depicted in Figure 7. The highest influencing parameter is density with 13.3% to 18.4%, followed by melting enthalpy with 10.8% to 13.6%, heat conductivity with 7.9% to 9.8% and melting temperature with 5.5% to 7.8%. All other material properties have always an influence of less than 6.0% and C I has almost no impact. For ρ, λ, L, c s and especially c l the influence is higher during discharging than during charging and for T m it is vice versa. Again, the impact is, in general, somewhat higher when lowering the parameters than when increasing them. The maximal absolute relative influence of a variation in material properties on the charging and discharging power over time is depicted in Figure 7. The highest influencing parameter is density with 13.3% to 18.4%, followed by melting enthalpy with 10.8% to 13.6%, heat conductivity with 7.9% to 9.8% and melting temperature with 5.5% to 7.8%. All other material properties have always an influence of less than 6.0% and has almost no impact. For , , , and especially c the influence is higher during discharging than during charging and for it is vice versa. Again, the impact is, in general, somewhat higher when lowering the parameters than when increasing them.   The maximal absolute relative influence on the global liquid fraction is presented in Figure 8 as a function of the variation in material properties. The highest influencing parameter is density with a variation of 4.3% to 4.8%, followed by melting enthalpy with 3.3% to 4.2%, melting temperature with 2.5% to 3.7%, and thermal conductivity with 2.6% to 3.6%. For the other parameters, the influence is always below 2.0% and C I , η and β as well as c l during charging and c s during discharging have almost no impact.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 12 of 25 The maximal absolute relative influence on the global liquid fraction is presented in Figure 8 as a function of the variation in material properties. The highest influencing parameter is density with a variation of 4.3% to 4.8%, followed by melting enthalpy with 3.3% to 4.2%, melting temperature with 2.5% to 3.7%, and thermal conductivity with 2.6% to 3.6%. For the other parameters, the influence is always below 2.0% and , and as well as during charging and during discharging have almost no impact. By analyzing the results as a whole, the material properties can be divided into three groups:

•
Material properties that have no, or only little effect, if they are varied like indicated here: only has an influence that is always below 1% (maximum of 0.12%). Here, one has to keep in mind that is varied by several orders of magnitude in literature, and not only by ±10%.

•
Material properties that have an influence on the power and the global liquid fraction, but none on the mean power (i.e., the total heat transferred): the material properties , , and − mainly influence the maximum deviations and not the mean power. The reason why also has a slight influence on the average output is that the storage system can still absorb or release a certain amount of sensible heat at the end of charging and discharging.

•
Material properties that also influence the mean power (i.e., the total heat transferred): all remaining material properties ( , , , , ) also affect the mean power when varied, but , and only slightly.
The power over time curve during charging is shown in Figure 9 for the reference case and an increased melting enthalpy. Throughout the whole simulation, the deviation is positive-the mean power is increased. In contrast, Figure 10 shows that an increase in thermal conductivity (also for charging) leads to a different type of deviation. First, the power is increased for the variation compared to the reference, but later it decreases-the mean power is only affected slightly. The same phenomenon can be seen for discharging and the contrary effect can be seen for lowering the melting enthalpy and heat conductivity ( Figures A1 and A2 in the Appendix A). This behavior is, in general, identical for all material properties that either affect only the maximum deviation of the power and the global liquid fraction or also the mean power. By analyzing the results as a whole, the material properties can be divided into three groups:

•
Material properties that have no, or only little effect, if they are varied like indicated here: only C I has an influence that is always below 1% (maximum of 0.12%). Here, one has to keep in mind that C I is varied by several orders of magnitude in literature, and not only by ±10%. • Material properties that have an influence on the power and the global liquid fraction, but none on the mean power (i.e., the total heat transferred): the material properties η, β, λ and T l − T s mainly influence the maximum deviations and not the mean power. The reason why λ also has a slight influence on the average output is that the storage system can still absorb or release a certain amount of sensible heat at the end of charging and discharging.

•
Material properties that also influence the mean power (i.e., the total heat transferred): all remaining material properties (ρ, L, c l , c s , T m ) also affect the mean power when varied, but c l , c s and T m only slightly.
The power over time curve during charging is shown in Figure 9 for the reference case and an increased melting enthalpy. Throughout the whole simulation, the deviation is positive-the mean power is increased. In contrast, Figure 10 shows that an increase in thermal conductivity (also for charging) leads to a different type of deviation. First, the power is increased for the variation compared to the reference, but later it decreases-the mean power is only affected slightly. The same phenomenon can be seen for discharging and the contrary effect can be seen for lowering the melting enthalpy and heat conductivity ( Figures A1 and A2 in the Appendix A). This behavior is, in general, identical for all material properties that either affect only the maximum deviation of the power and the global liquid fraction or also the mean power.

Boundary and Initial Conditions
Finally, the outcome of the variation of the initial and boundary conditions is presented and again an overview, including the sign of the resulting changes in the mean power is shown in the Appendix A in Table A3. When analyzing these results, one has to keep in mind that in general the parameters are varied by ±1 K and ±10%, but the heat losses, the inlet temperature (constant value vs. curve) and the inlet velocity (constant vs. parabolic profile) are varied in an on/off manner. Therefore, the influence is in general in a plus and minus direction, and only in one direction for the parameters with on/off variations. Consequently, the later presented ranking is based on a comparison of a mean of the variation by +1 K and −1 K or +10% and −10% with a variation from off (e.g., constant value) to on (e.g., profile). Figure 11 shows the mean deviation of the charging and discharging powers. The highest absolute influence can be seen for the heat losses with 21.4% and 16.4% during charging and discharging, respectively. The initial temperature and the inlet temperature have an influence of 0.8% to 2.2% on the mean power when varied by ±1 K. All other variations of the boundary conditions (inlet velocity, inlet temperature curve and inlet velocity profile) have almost no influence on the mean power.

Boundary and Initial Conditions
Finally, the outcome of the variation of the initial and boundary conditions is presented and again an overview, including the sign of the resulting changes in the mean power is shown in the Appendix A in Table A3. When analyzing these results, one has to keep in mind that in general the parameters are varied by ±1 K and ±10%, but the heat losses, the inlet temperature (constant value vs. curve) and the inlet velocity (constant vs. parabolic profile) are varied in an on/off manner. Therefore, the influence is in general in a plus and minus direction, and only in one direction for the parameters with on/off variations. Consequently, the later presented ranking is based on a comparison of a mean of the variation by +1 K and −1 K or +10% and −10% with a variation from off (e.g., constant value) to on (e.g., profile). Figure 11 shows the mean deviation of the charging and discharging powers. The highest absolute influence can be seen for the heat losses with 21.4% and 16.4% during charging and discharging, respectively. The initial temperature and the inlet temperature have an influence of 0.8% to 2.2% on the mean power when varied by ±1 K. All other variations of the boundary conditions (inlet velocity, inlet temperature curve and inlet velocity profile) have almost no influence on the mean power. Figure 11. Absolute relative variation of the mean power for charging and discharging when varying the boundary and initial conditions. In Figure 12, the maximum deviations of the charging and discharging power are presented for the variations of the initial and boundary conditions. The highest influence can be seen for the heat loss with 43.3% and 56.3% followed by the inlet velocity profile with 43.7% and 53.4% and the inlet velocity with 33.8% to 42.2%. The variation of the remaining initial and boundary conditions has an influence of 9.3% to 18.0%. The power during the charging and the discharging is depicted in the Appendix A in Figures A3 and A4 respectively for the reference case and the one with heat losses. As a comparison, the charging process with an inlet velocity variation of +10% is also examined in the Appendix A in Figure A5. In general, for all initial and boundary condition variations-except for the heat losses-the largest deviations occur only for a very short time interval at the beginning of the charging and discharging process (see e.g., Figure A5). This is also the reason why they have only a rather small effect on the average performance. Figure 11. Absolute relative variation of the mean power for charging and discharging when varying the boundary and initial conditions. In Figure 12, the maximum deviations of the charging and discharging power are presented for the variations of the initial and boundary conditions. The highest influence can be seen for the heat loss with 43.3% and 56.3% followed by the inlet velocity profile with 43.7% and 53.4% and the inlet velocity with 33.8% to 42.2%. The variation of the remaining initial and boundary conditions has an influence of 9.3% to 18.0%. The power during the charging and the discharging is depicted in the Appendix A in Figures A3 and A4 respectively for the reference case and the one with heat losses. As a comparison, the charging process with an inlet velocity variation of +10% is also examined in the Appendix A in Figure A5. In general, for all initial and boundary condition variations-except for the heat losses-the largest deviations occur only for a very short time interval at the beginning of the charging and discharging process (see e.g., Figure A5). This is also the reason why they have only a rather small effect on the average performance.
The maximal deviations of the global liquid fraction are shown in Figure 13 when varying the initial and boundary conditions. The heat loss has the greatest influence with 8.2% to 13.2%, followed by the inlet temperature with 2.4% to 2.8%. All other variations of the initial and boundary conditions have a minimal impact of 0.1% to 1.2% on the maximal deviation of the global liquid fraction. The maximal deviations of the global liquid fraction are shown in Figure 13 when varying the initial and boundary conditions. The heat loss has the greatest influence with 8.2% to 13.2%, followed by the inlet temperature with 2.4% to 2.8%. All other variations of the initial and boundary conditions have a minimal impact of 0.1% to 1.2% on the maximal deviation of the global liquid fraction.
Lastly, a comparison between experimental and numerical results with and without heat losses is performed on the mean power during the charging and discharging process until 7000 s. For charging, the achieved absolute mean powers are 26.2 W, 18.3 W, and 20.7 W for the experiment, the simulation without heat losses, and the simulation with heat losses, respectively. In the same order, the absolute mean powers during discharging are 14.1 W, 15.6 W, and 14.5 W. This leads to an absolute relative deviation of the simulation with heat losses from the experiment of 20.8% and 2.5% for charging and discharging, respectively. This is a distinctive decrease compared to the numerical results without heat losses; here the absolute relative difference was 29.9% for charging and 10.6% for discharging.   The maximal deviations of the global liquid fraction are shown in Figure 13 when varying the initial and boundary conditions. The heat loss has the greatest influence with 8.2% to 13.2%, followed by the inlet temperature with 2.4% to 2.8%. All other variations of the initial and boundary conditions have a minimal impact of 0.1% to 1.2% on the maximal deviation of the global liquid fraction.
Lastly, a comparison between experimental and numerical results with and without heat losses is performed on the mean power during the charging and discharging process until 7000 s. For charging, the achieved absolute mean powers are 26.2 W, 18.3 W, and 20.7 W for the experiment, the simulation without heat losses, and the simulation with heat losses, respectively. In the same order, the absolute mean powers during discharging are 14.1 W, 15.6 W, and 14.5 W. This leads to an absolute relative deviation of the simulation with heat losses from the experiment of 20.8% and 2.5% for charging and discharging, respectively. This is a distinctive decrease compared to the numerical results without heat losses; here the absolute relative difference was 29.9% for charging and 10.6% for discharging.  Lastly, a comparison between experimental and numerical results with and without heat losses is performed on the mean power during the charging and discharging process until 7000 s. For charging, the achieved absolute mean powers are 26.2 W, 18.3 W, and 20.7 W for the experiment, the simulation without heat losses, and the simulation with heat losses, respectively. In the same order, the absolute mean powers during discharging are 14.1 W, 15.6 W, and 14.5 W. This leads to an absolute relative deviation of the simulation with heat losses from the experiment of 20.8% and 2.5% for charging and discharging, respectively. This is a distinctive decrease compared to the numerical results without heat losses; here the absolute relative difference was 29.9% for charging and 10.6% for discharging.

Comparison with Previous Studies
In this section, the outcome of the current study is compared with results from literature. It has turned out that a simple one-to-one comparison should not be done, as the initial variables, the applications and the number of varied parameters as well as their variation range may be very different. In addition, Dolado et al. [27] and Mazo et al. [28] used a hypercube sampling and assumed a normal distribution of the varied parameter. In contrast, the analysis of Zsembinszki et al. [30] and the study presented here are based on a simple variation with respect to a reference case. Nevertheless, a comparison will still give more insights into the interdependencies and may reveal details for the need for future research. Therefore, in Table 5, the most influential parameters are listed together with their assumed variation range for all studies mentioned above. On average, T m is the most influential parameter in the studies from literature. Whereas, the heat loss is the most dominant parameter in this study. It is worth mentioning that the effect of heat loss is not investigated in the studies found in literature and might be distinct smaller due to proper insulation. A comparison of all studies that are investigating the influence on the mean power shows that the ranking of the most influential material properties is still different even though the influenced parameter is identical. In the present study, ρ and L are the most relevant material parameters and in the literature these are T m and L [27] or ∂h ∂T and ρ as well as T m [30] (the ranking for the deviation of the mean power is somewhat different from the ranking based on the deviation of the mean absolute power). Dolado et al. [27] also studied, among other parameters, the time until the output temperature reaches a certain temperature threshold. Here, the ranking is again different from the one when looking at the mean power. The ranking given by Mazo et al. [28] is different for the heating energy demand than for the cooling demand and the ranking depends on the PCM layer thickness, too. The corresponding ranking in Table 5 is the one given in the conclusion section of the original paper [28]. Table 5. Comparison of the most influential parameters and their assumed variation.

3. 4. 5.
Dolado et al. [27] mean power deviation Zsembinszki et al. [30] mean absolute power deviation Mazo et al. [28] saving in yearly energy consumption T m ±1 K * For the mean absolute power deviation in Zsembinszki et al. [28], ρ and ∂h ∂T are equally important; ** For the mean power deviation in Zsembinszki et al. [28], T in , ρ and ∂h ∂T are equally important; *** For the mean power deviation in Zsembinszki et al. [28], c HTF is equally important as ρ HTF .
Not only the ranking of the most important parameters varies between and within the different studies but also the overall impact of the variations varies widely. Mazo et al. [28] report that the yearly energy demand for cooling the studied test cell changes by about ±10% when all parameters are varied together. Moreover, for a variation of all parameters, Dolado et al. [27] have found that this causes differences in the mean power of about ±3.75%. In contrast, Zsembinszki et al. [30] state that changing T in by ±1 K already leads to a shift of the mean absolute power by about ±11%. Similarly, the results of the current study revealed that a variation of ρ by ±10% leads to changes in the mean power of about ±8%, and accounting for the heat losses or not leads to a difference in the mean power of 19.4% and 13.5% for charging and discharging, respectively.

Conclusions
Within this paper, a basic sensitivity analysis is performed to understand the sources of uncertainties within a CFD simulation of a vertical shell-and-tube LTESS during charging and discharging. First, the influence of several numerical parameters was investigated, and it was concluded that the major part of the uncertainties must be found elsewhere. Consequently, a sensitivity analysis of the material properties, as well as boundary and initial conditions and measurement positions, was carried out. In order to keep the simulation effort within limits, a simple analysis was performed in which the parameters were varied by a certain value in relation to a reference case. Ten material properties were varied by either ±10% or ±1 K, and six boundary and initial conditions by either ±10%, ±1 K, or in an on/off manner, and the resulting deviation of the mean power, as well as the maximum deviation of the power and the global liquid fraction, were analyzed. In addition, the position of the temperature measurements within the simulation was varied by ±1 mm. The most important results are: • The variation of the material properties had a very large effect on the power and the global liquid fraction. In conclusion, uncertainties in the material properties may affect the results of a CFD simulation of a LTESS significantly.

•
The material properties can be subdivided into three groups: Material properties that have no or only little effect: only C I has an influence that is always below 1% (maximum of 0.12%), but one has to keep in mind that it is varied by several orders of magnitude in literature. Material properties that have an influence on the power and the global liquid fraction slope, but none on the mean power: the material properties η, β, λ and T l − T s mainly influence the maximum deviations and not the mean power. Material properties that also influence the mean power: all remaining material properties (ρ, L, c l , c s , T m ) also affect the mean power when varied, but c l , c s and T m only slightly.
• . Q loss (up to 21.4%), ρ (up to 9.2%) and L (up to 6.4%) have the largest influence on the mean power.

•
The ranking of the parameters and the magnitude of the influence is sometimes distinctively different in studies found in literature.

•
Comparing the maximum deviation of the power can be misleading as the deviation may only occur as a short peak, and has no effect on the overall performance.

•
For the investigated storage configuration and orientation, a radial variation of the temperature measurement position had a much higher impact than a vertical variation. A displacement of the measuring point in radial direction by ±1 mm already led to a maximum deviation in the measured temperature of more than 5 K.
The above points indicate that the results of numerical models of latent thermal energy storage units are subject to substantial uncertainties. This is to a large extent a direct consequence of the uncertainties in the applied material properties and the imposed boundary and initial conditions (the order of magnitude of the variation agrees for many parameters with the uncertainties/variations found in literature [3,27,28,30]). In general, however, these uncertainties are not highlighted, but rather hidden. A fair comparison with experiments, including uncertainties like it is performed by Dolado et al. [27] and Mazo et al. [28], is hard to find in literature. One reason for this is certainly the large number of required simulations, which can be an exclusion criterion, especially for CFD studies.
A first next step should include detailed analyses of standard validation experiments. These analyses need to involve the actual uncertainty of the material properties, their temperature dependency, and heat losses, as well as the uncertainties of the experiment. To ensure that the whole parameter range is taken into account, sophisticated sensitivity analysis methods need to be applied. Moreover, the error caused by the Boussinesq approximation and the influence of the chosen switch-off method should be examined in detail. Funding: This publication was funded by the German Research Foundation (DFG) and the University of Bayreuth in the funding program Open Access Publishing.  Table A3. Overview of the influence on the magnitude of the mean power for all varied material properties as well as boundary and initial conditions. The charging and discharging power is always defined positive.  Figure A1. Power over time curve for charging of the standard case and the case with a decreased melting enthalpy as well as the regarding relative deviation. Figure A1. Power over time curve for charging of the standard case and the case with a decreased melting enthalpy as well as the regarding relative deviation. Figure A1. Power over time curve for charging of the standard case and the case with a decreased melting enthalpy as well as the regarding relative deviation.