**1. Introduction**

In Europe, the food and beverage industry is the fifth largest energy consumer [1], consuming the highest share of low-temperature heat demand of the total process heat demand [2]. Using renewable energies and increasing energy efficiency are ways to reduce greenhouse gas (GHG) emissions [3]. First, the potential to reduce process energy demand can be identified using process integration methods [4]. If the possibilities of increasing the energy efficiency at the process level through direct heat recovery are exhausted, further energy saving is possible by means of total site heat integration [5], with a specific focus on low-temperature processes [6] and heat recovery loops (HRL) [7].

In many industries, such as the dairy and beverage industry, the individual processes are operated non-continuously for product changes and cleaning reasons. This poses additional challenges to heat recovery, often requiring the use of energy storage. Two approaches to characterizing batch processes and its flows are the time average model (TAM) and the time slice model (TSM). The TAM indicates heat recovery potential, which can be achieved by indirect heat transfer and storage based on an HRL. Olsen et al. [8] present a systematic thermal energy storage integration method based on the TAM approach. Their indirect source and sink profile (ISSP) method enables the concrete identification of suitable sink and source profiles and the dimensioning of intermediate circuits and heat storage [9]. Due to the integration of multiple sink and source streams, the utilization of the storage capacity increases the overall availability and performance of sources and sinks. For sensible heat storage, there are closed stratified tanks (ST) and fixed temperature variable mass (FTVM) open storage tanks. FTVM tanks require a larger storage volume than ST tanks, but they are easier to operate and suitable for very large circulating volumes and small temperature differences. ST tanks need only half the volume and can operate at temperatures above 100 ◦C when placed under pressure. Due to the temperature dependency of fluid density, higher temperature fluid (lower density) forms a hot zone at the top of the tank, while lower temperature fluid forms one at the bottom. These zones are separated by a narrow transition area (between *h*Tu and *h*Tl)—the thermocline. For varying source and sink streams, the volume of fluid in each temperature zone increases or decreases, moving the thermocline vertically within the storage tank (see Figure 1). A completely loaded or unloaded storage tank, in addition to other mixing effects, results in a loss of stratification, loss of heat recovery, and even temporary shutdown of the HRL [10].

**Figure 1.** Heat recovery loop (HRL) based on a stratified tank (ST) and thermocline control.

Dynamic simulation is necessary for evaluating and proving the performance of an HRL system in any scenario. This is due to arbitrary production and logistic timing, the changing operation time of process states, thermal inertias, as well as dynamic interaction. Past models have focused on either the heat recovery side of the HRL [11] or the ST unit. These models assumed that there was a perfect interface between the cold and hot area of the ST [12], which is a simplification of practical installations. Walmsley et al. [13] studied an experimental scale model of an industrial ST, including the interaction between flows and thermocline movement and growth, but did not incorporate this element into their HRL models. The growth and movement of the thermocline accelerate the loss of stratification, reducing the effective heat recovery capacity of the tank over time and leading to reduced energy savings. Baeten et al. [14] presented a model that incorporates the buoyancy and mixing in one-dimensional stratified storage tanks in building energy simulations to estimate the uncertainty of parameters as a function of the storage capacity caused by the storage tank model. Campos Celador et al. [15] showed that the storage tank model used in such simulations can have major effects on calculated annual savings and design decisions. Powell and Edgar [16] introduced an adaptive-grid model with a high-resolution grid in the region of thermal stratification. This ensures higher accuracy and simulation performance than other one-dimensional models, and has speed advantages over computational fluid dynamics (CFD) models in repetitive simulations (e.g., for optimization, validation, control, prediction). Schlosser et al. [17] demonstrated a multi-node (MN) approach can represent a good compromise between accuracy and simulation time for the high-fidelity simulation of HRL systems. Numerical diffusion errors were reduced by approximately 35% using a variable layer height (VLH) model. The stratification of the physical thermocline growth could be re-established while still operating the ST by siphoning the mixed thermocline area from the tank [13]. Furthermore, a control system consisting of a single flow and system control for maintaining stratification was applied. The system control operated depending on the thermocline position, enabling restratification after being completely charged or discharged and a single flow control adjusted the individual streams to the target temperatures *Thot* or *Tcold*. This previous simulation study [17] included a typical week of production for tank sizing and evaluation, taking no account of stochastic changes (e.g., order quantity) or deterministic influences, such as recipe sequences. As an extension, the present work simulates several weeks of production, based on synthetic time series with probabilistic influence (e.g. changing recipes and order amounts) by Monte Carlo (MC) simulation to generate a sensitivity evaluation of the HRL system. Atkins et al. [10] highlighted the optimal sizing of an ST as a prerequisite for economic efficiency and operability. Moreover, Meschede et al. [18] applied MC methods to highlight the risk caused by using a single reference time series for designing energy systems. The MC method helps to define the mean of the heat recovery rate (HRR) for different tank sizes more precisely, and to better assess the risk of fluctuating heating and cooling load due to stochastic influences [19]. A prerequisite for this, however, is an accurate and high-performance VLH model.

The aim of the paper is to improve the design and operation of industrial HRLs by constructing a comprehensive model of storage and heat exchanger operations, simulated with a MC approach. Therefore, the first step of this study is to validate the performance and accuracy of the VLH model. Secondly, the effect of the variability of streams on the sizing approach in terms of heat integration potential and operability is evaluated by MC simulation, with real process data from the dairy industry. The MC method considers the sensitivity of stochastic influences (order quantity, sequence of recipes) and deterministic relations (production process: idling, operation, cleaning in place (CIP)) on the volatility of sink and source profiles. The results are probability distributions of the HRR for different storage dimensions of the present case study.

The article essentially contains two objectives: (1) the validation of the accuracy and performance of the VLH multi-node (MN) model and (2) a risk assessment of the robustness of the ST design and control by MC simulation, based on the validated model. The validation of the VLH model, developed by Schlosser et al. [17], is carried out on the basis of the validation experiments introduced by Walmsley [20] in comparison with a real laboratory tank and a CFD model (Sections 2.1 and 3.1). The simulation of a VLH model-based HRL, designed according to the ISSP method (Section 2.3), is realized by stochastically generated time series as input following the MC approach (Section 3.3). The generation based on cumulative distribution functions (CDF) was developed in the context of this work (Section 2.2) and applied (Section 3.2) for measurement data of a case study [12]. The results of Section 3 are summarized in Section 4.

#### **2. Materials and Methods**

The following sections describe the methods for validating the VLH model and the evaluation approach for the HRL dimensioning. The first section presents the different validation experiments proving the accuracy and performance of the VLH model. Subsequently, the generation of stochastic time series based on random influences and deterministic correlations is described. Finally, the evaluation approach of a robust storage sizing based on dynamic simulation is presented.

#### *2.1. Verification and Validation of the Variable Layer Height Model*

Verification and validation tests confirm if a model reproduces a certain (physical) behavior qualitatively and quantitatively correct, or with tolerable deviations. In this context, verification means checking whether the model correctly implements the derived equations. Verification can be understood as a qualitative examination of the transformation and the functionality of the model. Validation quantitatively checks the representation of the original system by the model. The quantitatively correct representation of the reference system is determined by subjectively defined tolerances, so that a model can be described as sufficiently accurate within a certain error deviation [21].

The model is examined in two dimensions–performance and accuracy. To calculate a large number of probable time series for a production week in an appropriate time frame, sufficient simulation speed is required. In terms of verification, the functionality is given by using universal physical descriptions. The validation is based on data from a CFD model, which is validated by a lab-scale, insulated ST (6.44 L) made of Perspex [20]. In this way, sufficient accuracy of the VLH model can be demonstrated by various experiments. After the numerical diffusion errors are minimized [17], the position of the thermocline, characterizing the heat recovery capacity, is used as a measure of the error deviation. As a validation quality, the following unitless key metrics are developed. The percentage of ideal cases (PIC) represents the share of the ideal ST case *A*0 compared to the areas of cold *A*1 and hot *A*2 losses of stratification, as illustrated in Figure 2. PIC quantifies the level of stratification. The PIC number is defined as:

$$PIC = 1 - \frac{A\_1 + A\_2}{\left(\frac{A\_0}{2}\right)}\tag{1}$$

**Figure 2.** Illustration of variables used to define the validation criteria in a ST.

This results in the following operating ranges [20]:


The position and height of the thermocline *htm* gives information about the available heat recovery capacity. The dimensionless height *htm*,*<sup>i</sup>* of the thermocline is determined by the division of the actual height *htm* and the total tank height *H*. Associated with the thermocline thickness *Ht* the degradation of the stratification is quantified.

$$H\_t = h\_{tu}(T\_{tu}) - h\_{tl}(T\_{tl}) \tag{2}$$

This work is now based on the validation parameters of the PIC, characterizing the degree of thermocline degradation. The deviation of the PIC from the experimental values in different experiments is quantitatively assessed. The height of thermocline *htm*,*<sup>i</sup>* is qualitatively discussed by means of a targeted loading and unloading profile. The following validation measures are conducted. The settings and boundary conditions of each experiment are presented in Table 1.



**Table 1.** Settings and boundary conditions of validation.

#### *2.2. Generation of Stochastic Time Series for Monte Carlo Simulation*

Understanding the variability of the available sources and sinks is prerequisite for an optimized volume of the required thermal storage and to achieve a balanced system. As already outlined, several different factors have an influence on the height and temporal course of the heat source and heat sink processes. Within the framework of a simulation study, reliable and comparable energy savings considering the dynamic system's behavior, can be calculated. A sensitivity analysis considers the influence of fluctuating factors on the target variables robustness of the HRL dimensioning and HRR. If probabilities are assigned to the stochastic influencing factors, time series that occur with a certain probability can be generated from the deterministic correlations. The generated time series serve as input data for the dynamic simulation. Thus, the simulation validates the resulting heat recovery potential, depending on the tank size, regarding probability. The methodology for generating time series used in this paper follows the steps shown in Figure 3.

**Figure 3.** Flow chart of the approach for generating weekly time series of the heat source and sink processes.

The first step is to collect data. These data are usually measured energy flows and production quantities. Data from a case study in New Zealand are used for this purpose [12]. The analysis of all essential influencing parameters of a production week requires a recording of the measured data with a suitable high temporal resolution. Seasonal effects (e.g., the weather, or availability of milk) is neglected for the heat recovery problem. If no measurement data is available (e.g., for greenfield planning), the necessary parameters could be assumed as normally distributed around average values from literature data or expert opinion.

In the second step, the time series are analyzed for possible deterministic correlations. Since the flow rate of the streams varies due to process fluctuations and discontinuous operation of the plants, which are caused by regular cleaning, maintenance work, plant failure, product change, and changing milk supply. Supply temperatures are generally quite constant, due to robust process control, and do not represent a major cause of variability in this case. Thus, the sequence of idling, production, and CIP is mandatory.

Moreover, the probabilistic components of the weekly values are described by CDF based on cumulative frequency plots. Based on the result of the statistical analysis, multiple weekly time series are generated in step three.

Due to the interaction of product changes and cleaning constraints on the energy demand, the streams have different demand profiles. Generally, two operating states can be distinguished, especially in the food industry: "on-production", associated with an energy demand, and "off-production" (CIP or breakdown), without any energy demand. Duration, frequency, and amplitude (i.e., production rate) of these operating conditions can be determined by statistical analysis of process time series. Constant and fluctuating energy demands can be detected. Based on that, the following three characteristic time series can be derived.

Using dairy case study data [11], the method for generating time series is demonstrated. The measurement data of two processing months are statistically evaluated to derive the weekly demand time series. Based on the load profiles in Figure 4, the frequency distributions of the heat flow rate and of the durations *t* for the operating states "on-production" and "off-production" are represented in the form of histograms (Figure 5a). The determinations of the frequency and cumulative probability distributions are shown in Figure 5, using the time series type (b) for stream Whey A, as an example. The computations are carried out with MATLAB®(9.3.0.713579 (R2017b), The MathWorks Inc., Natick, MA, USA).

**Figure 4.** Characteristic time series for a (**a**) constant and (**b**) fluctuating heat flow with time intervals in off-production and (**c**) a continuously fluctuating heat flow.

**Figure 5.** (**a**) Histograms and (**b**) empirical cumulative distribution functions (CDF) of the heat flow rate, duration of "on-production" and duration of "off-production" for the example of Whey A.

The generation of the time series begins with the decision of whether production is "on" or "off". The probability for this decision is calculated by the mean value from the mean duration of the "on" and "off" periods. The CDF based on the cumulative frequency plots of the heat flow rate (Figure 5b) over the production period is used to set probable heat flow rate for every time step during the operating state "on-production", by which the demand is scattered with the corresponding variance. The distributed durations of the operating states *ton* for "on-production" and *toff* for "off-production" can also be determined from the CDF functions (Figure 5b).

In this context, MC simulations allow the consideration of numerous probabilistic inputs to determine the distribution of probabilistic results. In line with previous studies [22–24] on the robustness of energy systems, probabilistic times series based on probability density functions (PDF) and CDF are used. In this way, representative data can be generated for multiple purposes, e.g., validation and sensitivity analysis.

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

The stochastic time series of the heat sources and sink processes are used as input data for the thermal simulation of the HRL from Figure 1. In the sensitive simulation study, the tank size is designed for the average stream data according to the ISSP method. Furthermore, the target temperatures *Thot* and *Tcold* of the intermediate circuits will be varied. Finally, the HRRs for different scenarios are presented depending on their probability of occurrence. The HRL (Figure 1) simulation model based on MATLAB/Simulink® consists of heat exchangers, a stratified storage tank, and a control system. The heat exchangers are modelled according to the ε-NTU [25]. The VLH model based on the MN approach is chosen as the modelling approach for the ST [17]. The system control is divided into two parts, as shown below.

The position and width of the thermocline describe the hot and cold zone capacity and therefore serve as control criteria in the system control. The system control (Figure 6a) specifies various operating states (i.e., normal operation, loading, unloading) depending on thermocline location in the ST. If the middle of the thermocline reaches the lower or upper end of the storage tank, this is considered to be completely loaded or unloaded. The source and sink circuits are then deactivated to regenerate the stratification. For hysteresis reasons, the sources and sinks remain deactivated until the corresponding zone again occupies 10% of the storage volume. The single mass flow control (Figure 6b) adjusts the streams to the target temperatures *Thot* or *Tcold* by mass flow control. Moreover, the growth of thermocline is contained to a certain extent *Ht*,*max* by siphoning while the HRL and ST are still in full operation.

**Figure 6.** (**a**) System controls depending on the thermocline position and (**b**) single mass flow control adjusting *Thot*and *Tcold*.

For evaluating the quality of the HRL system, an energetic target value for the maximum heat recovery (MHR) potential is determined according to the TAM method. The design of the HRL is carried out according to the ISSP method. Based on average stream data of a repetitive stream-wise repeat operation period (SROP) of one week, the hot and cold streams are defined. The stream data consist of average load, supply, and target temperatures. Considering a minimum temperature difference Δ*T*min, the hot and cold composite curves (CC) are formed. The overlapping area indicates the heat recovery potential (Figure 7a).

**Figure 7.** From (**a**) the indirect source and sink profile (ISSP) to (**b**) the HRL dimensioning.

Depending on the form and position of the CCs, the temperatures of the intermediate circuits, as well as *Thot* and *Tcold*, can be read directly, in addition to the number of intermediate circuits (Figure 7b). With known temperature difference (*Thot* − *Tcold*) between the hot and the cold side of the

tank and the maximum product of the difference between heat sources . *QSources* and heat sinks . *QSinks* for a certain time *t*, the water-based storage volume can be calculated.

$$V = \frac{\max\left(\left(\sum \dot{Q}\_{\text{Sources}} - \sum \dot{Q}\_{\text{Sides}}\right) \cdot t\right)}{1.16 \frac{kWh}{m^3 \cdot K} \cdot \left(T\_{hot} - T\_{cold}\right)} \tag{3}$$

The resulting tank size is highly dependent on the time of the occurrence of sinks and sources. Therefore, the design on the basis of one measuring profile can lead to a decrease of the HRR, or worse, to a standstill. A reliable estimation of energy savings and assessment of the tank size is only possible with a sensitivity analysis. The simulation of multiple weekly time series with random thermocline positions at the start of production evaluates the reliability and robustness of the chosen tank size. The model is also simulated applying different tank sizes. Finally, the probability distribution of the heat recovery potential for different tank sizes is evaluated based on the MC analysis.
