**3. Results**

In the following sections, the results of the validation, the HRL sizing, and the sensitive simulation study applied to the dairy case study are presented.

#### *3.1. Validation of the Variable Layer Height Model*

The experiments described in Section 2.1 are performed for the laboratory tank, the CFD, and the VLH multi-node (MN) model. For each experiment and each model, temperature curves are compared for characteristic points over the dimensionless storage height. The degree of destratification is represented by the PIC value.

(1) (a) **Static Mode:** Numerical diffusion effects are reduced to an acceptable degree, as shown in the previous study [13]. For the lab tests, the heat losses cannot be eliminated completely by the thermal insulation. For all experiments (Figure 11a–c), the thermocline grows with the simulation duration due to numerical diffusion effects. The deviation from the ideal stratification (black line) measured by the PIC values thus also increases, as shown in Figure 11d. Using the VLH model, the influence of numerical diffusion is even lower than when using the CFD model.

(b) Incoming or outgoing mass flows are still deactivated. Heat losses, on the other hand, are taken into account, which leads to the following destratification for different simulation times in Figure 8a–c. The destratification of the three models takes place to a similar extent in Figure 8d. Since no vertical heat conduction effects are taken into account for the MN approach, the stratification can only degrade evenly over the height (Figure 8b).


**Figure 8.** Destratification of the (**a**) lab, (**b**) VLH model and (**c**) CFD model due to heat losses for the **Static Mode**, compared based on the (**d**) PIC value.

**Figure 9.** Thermocline **Movement** for a targeted tank charging sequence in comparison to the ideal position for the (**a**) lab, (**b**) VLH model and (**c**) CFD model, compared based on the (**d**) PIC value.

Moreover, the scalability from lab-scale sizes to industrial-scale has to be shown. Therefore, the lab-scale tank and the VLH model have the same aspect ratio (height/width) as common industrial applications. Moreover, the stratification is influenced by the flow characteristics, especially the residence time in the tank. The maximum flow velocity for the laboratory tank was 0.002 m/s, which can also be found in industrial applications. The modelled ST should be in between this range. For this reason, sufficient scalability is given.

**Figure 10.** Influence of turbulent flows on stratification of the (**a**) lab, (**b**) VLH model and (**c**) CFD model during **Charging**, compared based on the (**d**) PIC value.

**Figure 11.** Comparison of the (**a**) lab tank, (**b**) variable layer height (VLH) multi-node (MN) model and (**c**) computational fluid dynamics (CFD)-modelled tank in terms of the (**d**) PIC value for an adiabatic **Static Mode**.

#### *3.2. Stochastic Time Series of Dairy Plant*

For the case study, representative process data of the milk processing industry have been assembled into an exemplary complete system. This results in load profiles of the sources and sinks for a mean SROP, as seen in Table 2.


**Table 2.** Stream table of dairy case study, producing milk powder and cream [11].

Using the approach for generating time series, explained in Section 2.2, various load profiles can be generated for a SROP, as presented in Figure 12 for the example of Whey A.

**Figure 12.** (**a**) Original time series of Whey A compared to the (**b**) generated one.

Profiles can be generated in the same way for all streams from Table 2. Figure 13 shows four examples of the accumulation of total sources and sink streams. This results in multiple total sink and source streams, which represent randomly arranged operating states and loads, as seen in Figure 13. These can overlap in any imaginable way.

**Figure 13.** (**a**) Original and (**b**–**d**) exemplary total available source and sink streams of one stream-wise repeat operation periods (SROPs).

#### *3.3. Evaluation of the Robust Heat Recovery Loop-Dimensioning via Monte Carlo Simulation*

Figure 14 shows the exemplary dimensioning of the HRL (temperature and tank size) by ISSP for the mean values. A tank volume of 1000 m<sup>3</sup> is calculated. Since the concurrency of the sources and sinks is highly variable, the sizes 50, 100, 300, 500, and 2000 m<sup>3</sup> are investigated by the sensitivity analysis. For all simulations, a height/diameter aspect ratio of 3/1 is chosen.

**Figure 14.** (**a**) Heat recovery potential according to the time average model (TAM) and (**b**) HRL sizing by the ISSP method.

The simulation tool uses the storage temperatures, size, and loop temperatures targeted from the ISSP design. Next, the HRL is simulated for different sink and source profiles dictated by the control system. The simulation of multiple time series for each tank size with random thermocline positions at production start takes place. The distribution of the resulting HRRs for the different storage sizes and simulation runs are shown in Figure 15. The HRR is calculated by the ratio of recovered heat to the target value. For the case study with a minimum temperature difference of 5 K, a heat recovery target of 1250 MWh corresponds to the overlapping area of the hot and cold CC in Figure 14a.

**Figure 15.** Distribution (bars) and CDF (stairs) of the heart recovery rate (HRR) for (**<sup>a</sup>**–**f**) different tank sizes from 50 m<sup>3</sup> to 2000 m3.

Regarding the number of runs that must be carried out to achieve a distinctive probability distribution, some approaches use a convergence criterion to determine the number of simulations. In this paper, 200 simulations are run to generate probable time series, since this number shows reasonable distributions within an appropriate computing time. Storage tanks with > 1000 m<sup>3</sup> (Figure 15f) do not achieve any further heat recovery. A high HRR can be achieved with much smaller tanks. Only the standard deviation increases significantly with decreasing size.

In total, all tank sizes in combination with the system control are suitable to achieve HHR over 90%. With decreasing storage volume at the same aspect ratio, tank height and diameter are reduced. As a result, the incoming and outgoing flow velocities temporarily rise above the limits of the validity range, which is avoided by limiting the mass flow by the control system. For this reason, the yields decrease further for smaller tank sizes.

## **4. Discussion**

Previous studies [17] showed that a VLH MN modelling approach reduces thermocline growth, due to numerical diffusion effects while performing with suitable execution time and accuracy. The conducted validation tests confirm this and show lower thermocline growth than a CFD model. However, it should be noted that vertical heat conduction effects are not considered in the MN approach. Furthermore, the validity range of the MN approach must be restricted regarding the flow velocity, since turbulent flow effects are not taken into account. However, these operating points are to be avoided anyway for a stable HRL operation. In addition, scalability from lab-scale sizes to industrial-scale can be assumed based on the same aspect ratio (height/width).

The recording of thermal energy data for the design and evaluation of HRLs often requires a good deal of effort and is not representative of the whole operating range. One way to compensate for this is to generate one's own time series data, also called data farming. The representativeness of the generated time series increases with the increase in basis data. Figure 5a shows a left-shifted distribution of the loads. The distributions of the durations for "on-production" and "off-production" states, on the other hand, scatter wider. If no measurement data is available, input data can be randomly distributed around expert values.

Recent study [19] had shown the relevance of risk evaluation due to the variability of streams by MC simulation methods when assessing the sensitivity of the results. The results of this study confirm this. Although the design methods provide a storage capacity of 1000 m3, the results for 50 m<sup>3</sup> only show a HRR reduced by less than 1%. Even the insignificant further variance of the values suggests a stable operation of the HRL in combination with the existing control system over a wide operating range. A prerequisite for the extensive simulation runs is an accurate and high-performance VLH model. Based on this work, an evaluation method for robust system sizing, stable operation control, and reliable energy-saving prognosis is developed, which is applicable also to other thermal systems. Future work will investigate the robustness of heat recovery systems based on heat pumps with STs.

#### **5. Conclusions and Outlook**

The validation shows suitable accuracy and performance of the VLH model with reduced numerical diffusion effects, restricted by the validity range for non-turbulent piston flow. Furthermore, a method is developed that generates stochastic time series for a sensitivity analysis with few input values. A suitable form of sensitivity analysis is the MC simulation assessing the risk. The robustness of the HRL system is aided by the control system. The method also avoids oversizing in any of the probable operation states. The simulations show that the tank size has only a small influence on the HRR, since the mean value varies only between 95.48% and 95.90% for tank sizes of 50–2000 m3. The standard deviations are very small (< 0.54%) but increase with decreasing storage size. This leads to the conclusion that a 50 m<sup>3</sup> tank reduces the HRR by less than 1%, with a minimally larger variance, compared to the 1000 m<sup>3</sup> and 2000 m<sup>3</sup> tanks. Furthermore, it is now possible to better forecast reliable energy saving.

The usage of the VLH model for a model predictive control strategy is planned for future works. In addition, an evaluation benchmark for the verification of a sufficiently large sample of time series for a convergen<sup>t</sup> result should be introduced.

**Author Contributions:** Responsible for the conceptualization, F.S.; methodology, F.S. and H.M.; validation, F.S., M.R.W.W. and M.J.A.; investigation, F.S.; data curation, F.S., T.G.W. and M.R.W.W.; writing—original draft preparation, F.S.; writing—review and editing, M.P., H.M., R.-H.P. and T.G.W.; visualization, F.S.; supervision, T.G.W.

**Funding:** This research has been supported by the EU Project "Sustainable Process Integration Laboratory—SPIL", Project No. CZ.02.1.01/0.0/0.0/15\_003/0000456 funded by EU "CZ Operational Programme Research, Development and Education", Priority 1: Strengthening capacity for quality research, in collaboration with University of Kassel (DE), Technische Hochschule Ingolstadt - Institute of new Energy Systems (DE), and University of Waikato (NZ).

**Conflicts of Interest:** The authors declare no conflict of interest.
