Modeling, Simulation and Optimal Operation of Multi-Extraction Packed-Bed Thermal Storage Systems †

: Solar thermal power technologies require storage systems to mitigate the natural variability of solar irradiation. Packed bed thermal storage systems (PBTES) offer a cost-effective solution using air as heat transfer ﬂuid and rocks as a storage medium. Compared to its alternatives, however, PBTES presents a limited ﬂexibility of operation due to the conventional unidirectional ﬂow, which involves the progressive reduction of the outlet temperature during discharge and thus lowers the thermodynamic efﬁciency of the power cycle. The present study summarizes the progress on the design and optimal operation of a novel multi-extraction PBTES, a project that aims at mitigating its typically poor operational ﬂexibility for solar power applications. To this end, a one-dimensional model with a high spatial resolution of a PBTES was developed, which includes four intermediate outlet points along the axial direction to investigate the beneﬁts of optimal extraction operation. In order to reduce the computational burden, a coarser model of the storage system is used in combination with non-linear model predictive control (NLMPC). Through the optimal manipulation of the extraction valves, the output temperature is maintained close to a prescribed temperature throughout the discharge. The control admits not only constant temperature targets, but also time-varying scheduled proﬁles. This work describes the limitation of such a design and control approach and sets the direction for the future, more detailed analyses needed to demonstrate its applicability.


Introduction
Packed bed systems are considered a promising energy storage option to overcome the variability of renewable energy sources in applications like solar-thermal power generation [1][2][3][4][5] or wind power [6]. According to these different authors, the main advantages over competing technologies are higher operating temperatures, the possibility to use ambient pressure air as heat transfer fluid and to be placed after the receiver without additional heat exchangers, its good scalability, operational robustness and simplicity, as well as the low cost of the materials used.
Regarding power generation, the main shortcoming of this technology is the progressive degradation (decrease) of the heat transfer fluid (HTF) outlet temperature during the discharge process [4,7]. So far, the proposed technical solutions to stabilize this temperature during the heat retrieval have focused on the use of phase change materials (PCM), either used in combination with molten salts [8], or with packed bedrock systems [4,9]. Promising results regarding exergy efficiency and systems costs have been shown for the second option, especially compared with the complexity of cascade PCM systems [7].
Notwithstanding, the simplest solution to stabilize the outlet discharge temperature from a packed bed thermal storage systems (PBTES), other than increasing its size (mass), is by means of the selective use of the different thermal regions of the storage system. Each region is assumed to have a different temperature at the end of the charge, which could be used to produce a mixing flow with a desired temperature. This approach has been explored in the literature, as it will be detailed in Section 2. However, the importance of whether there exists an optimal strategy for the use of the thermal regions has been overlooked.
In PBTES, this selective use results in the extraction of mass flowrates that partially or totally affect the trajectory of the airflow within the storage volume. Extractions in radial or axial direction affect pressure conditions within the porous medium, which, in turn, modifies the flow path and the heat transferred throughout the system. While the heat transfer mechanisms that take place within current designs of bedrock storage systems are well understood, and there is little uncertainty in the response of the system given the usual uniaxial character of a constant airflow, this may not be the case when the airflow is intentionally modified in terms of direction and quantity. This is the starting point of a preliminary assessment of a novel control scheme that uses optimization tools to manipulate the flow of a multi-extraction packed bed storage system. It will be shown how, with such a configuration, it is possible to control the outlet temperature and thus increase the flexibility and efficiency of the system. The remainder of this paper begins with a review of the literature of packed bed systems (Section 2), followed by a section dedicated to materials and methods: design, modeling, control, and software used (Section 3). The model validation is subsequently described (Section 4), followed by several case studies and an extensive discussion (Section 5). The paper closes with some conclusions and an outline of future work (Section 6).

Literature Review
Generally, the design of PBTES consists of a cylindrical or conical bed of rocks with air flowing in opposite vertical directions when charging and discharging the system. Designs with horizontal uniaxial flow with trapezoidal sections have been also considered [5,6]. The main weakness of the uniaxial approach is the difficulty to maintain a constant temperature throughout long discharge periods. This shortcoming has been addressed in the literature by a combination of the following: increasing the size or thermal inertia, adding one or more layers of PCM, and segmenting or dividing the thermal storage system.
Michels and Pitz-Paal [10] demonstrated that arranging different PCMs in series (cascading), with rising fusion temperatures, increases the efficiency of the storage and maintains the desired output temperature for longer periods. A mixed approach using sensible and latent heat storage, a concept that has been demonstrated [9], has similar effects but at a lower cost.
Multi-storage or segmented systems, with the units connected in parallel or series (cascade), have been proposed in the literature to maintain a constant output temperature of thermal energy storage (TES) systems. These include sensible (rock, aluminum, etc.) and latent (PCM) cascade storage. For instance, Crandall and Thacher [11] proposed storage segmentation in order to avoid thermocline degradation when the temperature of the charging fluid (hot air from a solar field) decreases in 2004. Re-routing air at decreasing temperature to certain strata improves the efficiency of the charging process.
Segmentation was also mentioned in Kulakowski and Schmidt [12] in 1980 but with a different purpose, i.e., the use of parallel storage systems to obtain a more constant output temperature. They defined an algorithm to optimize the outlet temperature of a heat storage system including a by-pass line and an auxiliary heating system. The aim was to meet a consumer demand for hot air while maintaining a certain flow at a particular temperature. The optimal operation problem seeks to minimize the associated operational costs (back-up fuel). Although not implemented in their work, Kulakowski and Schmidt suggested that the same formulation can be applied to multi-storage systems, which would be connected in parallel, with respect to the fluid flow.
McTigue and White [13] reported on the benefits that segmentation could provide to thermal storage systems in the context of pumped TES, which included higher overall efficiencies, allowing for smaller particle diameters. More recently, these authors presented a comparative, including packed beds with radial and axial flow [14], where they concluded that radial flow schemes present similar efficiencies to non-segmented axial ones.
Regarding segmented packed beds operational benefits, Bindra et al. [15] proposed the so-called "sliding flow method" to improve the efficiency of segmented PBTES. The method has been used also in [16]. It consisted of charging the segments one by one until all of them are charged at the maximum charging temperature, thus avoiding stratification and reducing pressure losses. Their approach did not account for possible constraints in charging time to refill every segment, nor trade-offs regarding discharge temperature and discharge time. None of these aspects were taken into account in [11][12][13]15].
In light of the literature examined, which include significant efforts in the fields of modeling, simulation, an optimized plant design, it has been found that there exists room for improvement in the area of optimal operation of PBTES. Especially so with regards to segmented systems, where the possible combination of fluid flows at different temperatures to achieve a desired outlet temperature can be posed as an optimal constrained control problem. The present work completes previous efforts by introducing a novel theoretical design and operational framework that aims at improving the flexibility of packed bed storage systems. More precisely, this work answers the following research question: is it possible to design a simple model-based control algorithm to optimally select the flow path of a segmented PBTES so that a prescribed outlet temperature reference can be tracked? Figure 1 shows a schematic view of a packed bed storage system, showing the flow direction with variable section geometry, as well as the proposed multi-extraction solution. For a more detailed graphical representation of the non-segmented and segmented models, the reader is referred to [17] and [11] respectively. The proposed solution for bedrock TES system involves successive flow extractions in a controlled fashion, so that cost-effective operating strategies subject to constraints, e.g., the temperature of the mix or total mass flow rate, can be considered. It is assumed in this work, following the approach used in [15], the case in which the segments are physically detached so that a uniaxial flow can be considered within each segment.

Model Assumptions
Modeling and simulating fluid flows in porous media is generally a computationally expensive numerical problem. Following the example of other authors, a series of assumptions to reduce its complexity has been considered: variables depend on time and space in the axial coordinate; wall effects are neglected since D/dp > 40, following [4], after [18]; the intraparticle temperature gradient effect is non-negligible (Bi < 0.1), but its effect is simplified following [3], after [19], using a corrected number of transfer units (NTU) value; the inertial term (thermal capacity) in the energy balance for the fluid phase is neglected, following [20]; particles are considered in mono-dispersed and in homogeneous distributions [4]; no deformation or leakage [4] is considered; air thermal conductivity is independent from pressure, following [4], after [21]. The partial differential equation (PDE) system developed by Schuman [22] and modified by other authors yields, after simplification, a system of differential algebraic equations (DAE), as described in [4].

Packed Bed Model
Packed bed models generally rely on the spatial discretization of the porous medium. This paper follows the conventional approach of dividing the PBTES in layers or slices upon which balance equations are applied, and applies a similar numerical methodology as described in detail in [17]. Figure 2a illustrates an example of a network model with two partition elements between consecutive extractions, and five of these. Valves are represented by variable resistance elements R 1 , R 2 · · · R 5 , while the connecting branches are represented as R a , R b · · · R e . When the system operates in uniaxial mode, all valves except for R 5 remain closed. Figure 2b depicts an element, the variables assigned to it, and the relationships with the immediate downstream element. The validity of the simulated results depends significantly on the resolution of the resulting mesh [23]. The bedrock considered in this work is discretized in a spatial grid of 400 elements, with 80 layers or elements between consecutive extractions.
Equation (1) describes the relationship between the pressure drop at each element i as a function of the mass flow rate (trough the mass flux G =ṁ/A cross ) and the properties of the fluid, as well as the properties of the porous medium [3,24]. Energy balances for the fluid and solid volumes can be written, following [23], as: where h and e are, respectively, the enthalpy of the fluid and energy and solid media. As explained in [23], the time derivatives in Equations (2) and (4) can be solved using a forward Euler method, while the spatial derivative with a first-order forward difference.
Equations (2) and (4) are solved explicitly at each integration step. For each element i, the heat transferred to the fluid is proportional to the volumetric heat transfer coefficient α v , its volume, and the difference between the fluid temperature and the rock temperature of the element. The correlation of the volumetric heat exchange coefficient in Equation (6) is borrowed from [25].
Equations (2)-(5) are used to calculate the energy balance at each packed bed element for the fluid and rock phases. However, additional balance equations are needed to solve the networked model illustrated in Figure 2a. These include momentum equations for the valves and ductwork for the extractions and mixing sections, as well as energy balances to calculate the temperature of the mix, initial conditions and boundary conditions.
T(x = 0) = T in (11) p(x = 0, L) = p 0 , p L ∆ṁ(x = x 1 , x 2 , ..., Equations (7)-(9) express, respectively, mass-energy and momentum conservation, for every node i and branch j across the network. For every node, mixing at constant pressure in Equation (8) is considered. A simplified linear momentum equation along each branch of the system with negligible time derivatives is assumed in Equation (9). A set of initial and boundary conditions is required to solve the previous DAE. Equation (10) defines the initial conditions for the rock temperature, which will change depending on the scenario considered. The boundary conditions are summarized in Equations (11)- (13). Equation (11) refers to the inlet temperature, where the inlet section (x = 0) changes with the coordinate system, while Equations (12) and (13) define the pressure and massflow boundary conditions. As it will be shown in the next section, extraction flowsṁ ext,i may be imposed, which will result in an over-constrained problem. However, this can be relaxed by introducing the variables corresponding to the position of the valve illustrated in Figure 2a by the variable resistance elements.
Changes in fluid density affect the total pressure drop in the system. During charge for example, density decreases as the fluid temperature rises over time and, according to Equation (1), pressure drop increases. In order to reduce complexity, Equations (7)-(9) were solved assuming an average density of 0.6 kg/m 3 , which simplifies the assumption of constant fan pressure so that the prescribed flowrate can be obtained. The specific heat for air was assumed constant (C p,f =1050 J/kg K), while for the rock a linear interpolation of the data from [9] was used (C p,r = 800 + (1100 − 800)/(650 − 150)(T r − 150) J/kg K). Finally, the effective thermal conductivity was approximated by k eff = 0.5 + (2.0 − 0.5)/(650 − 150)(T r − 150) W/m K, based on results presented in [23].
The resulting differential-algebraic equation (DAE) system was solved using Simulink with the solver ode15s. For convenience, the origin of coordinates is flipped from the charge to the discharge stage in order to use the same code, which is a convenient approach also used in [17].

Optimal Operation
The simulation process can be divided into three steps: solving the flow and pressure system of equations, then the heat and temperature system of equations, and finally the integration of the rock temperature field. Figure 3 describes the flow of information used to simulate and optimally control the plant. In addition, the diagram shows the NLMPC block that calculates the optimal mass-flow extractions over the prediction horizon. The control problem can be described as minimizing the outlet temperature error subject to energy and mass balance constraints during the rolling horizon: The weights γ q can be selected to assign relative preference, for example when a price forecast is available for the electricity sold to the market or to a specific consumer. However, in the present work, the weights are selected in order to include (γ q = 1) or discard (γ q = 0) a future time interval of the receding horizon. This formulation is based on a simplified version of the epsilon-NTU method, where characterizes the effectiveness of the heat transfer process. It can be safely assumed that it approaches 1 in packed beds such as the one studied here. The prediction model divides the storage into five sections (n ext = 5) in the axial direction. Every iteration, the states are updated with the actual rock temperature, measured at specific points, to mitigate the effects of model inaccuracy.
Here, the temperature of the first 40 rock elements before each extraction is averaged for each section. The optimization problem defined by Equation (14) is designed to apply exclusively during the discharge phase, although it can be extended to cover also the charge phase.
The NLMPC optimization problem is solved with the active-set algorithm using the function fmincon (Matlab R , Mathworks, Natick, MA, USA) for non-linear, constrained optimization.

Model Validation
The validation of the model described in Section 3.2 is achieved by reproducing the experimental results of a reference system presented in [23,26], a concentrated solar power plant (CSP) located in Ait Baha, Morocco. The plant has a PBTES with 7.2 GWh, and a nominal electrical power output of 26 MW. Its design and operating parameters are summarized in Table 1.  Figure 4 shows the spatial temperature evolution at the end of cycles 1 and 30, starting from initial conditions of homogeneous temperature (150 • C). The duration of each cycle is 24 h, of which 8 h are used to charge the PBTES using the captured solar radiation. The rock temperature profiles are plotted at the end of the charge and discharge periods, replicating with high fidelity the results from [26], which demonstrates the ability of the model to reproduce the behavior of the reference plant under uniaxial flow conditions.  In the presence of extractions, the assumption of small temperature gradients in radial direction may not hold, and the model's ability to represent the real system would be questionable. However, if the cross-section is sufficiently small compared with the total length of the storage (D H), or if the extractions are effected in absence of radial flow (e.g., assuming physically detached segments), the approach presented in this work is expected to be valid. This is possible also because of the properties of model predictive controllers, which naturally mitigate the associated model inaccuracies by means of the feedback of state measurements [27].

Case Study
After validation in uniaxial flow conditions, several experiments are conducted to determine whether the proposed scheme achieves the desired outcome. This section is divided into two parts. First, a smaller, multi-extraction plant design is described and its results are compared with a larger reference packed bed, both under conventional uniaxial operation. Second, the extractions of the new design are optimized and the benefits discussed. The changes in the original design for the proposed plant assume that the ratio height-to-diameter has been multiplied by r 2 , where r is the volume ratio V/V 0 that takes the value 0.1, i.e., the new design is ten times smaller. The new dimensions are, therefore: D top = 27.4 m, D bottom = 21.9 m, and H = 5.4 m. For the numerical formulation to be valid, this geometry requires the assumption that the rock segments are physically divided in such a way that the extractions do not introduce radial flow. Figure 5 shows the respective output temperature in conventional operation mode, i.e., uniaxial mode. The resulting profiles evidence the shortcomings of using smaller storage volumes when no additional design or smart operating strategies are used, such as PCM or segmentation. Figure  5 also suggests that it is not possible to stabilize the outlet temperature through multi-extraction flows at the same levels as the original storage design, which averages around 620-630 • C. Given the low outlet temperature at the end of the discharge for the smaller storage design (below 400 • C), the maximum sustained outlet temperature by extraction mixing is likely to lay between 590 • C and 610 • C. This would ensure better partial load efficiency of the power block and lower moisture content at the outlet of the steam turbine [4] during the discharge. In order to verify these estimations, several experiments were carried out using the formulation introduced in Section 3.3. The control horizon, equal to the prediction (H c = H p ), decreases over the discharge period (16, 15, ..., 1), for this is where the NLMPC takes control over the operation. The values of the fixed hydraulic resistances are R a...e = [0.01, 0.01, 0.01, 0, 0] Pa/kg s −1 (Figure 2a). The value of the control resistances R 1...5 are computed every iteration once all extraction flows are known (Equations (1), (7), (9)-(13)). The length of each interval, ∆t, is set to one hour. The TES inlet temperature (T inlet,dis ) is assumed constant, which may not be appropriate under realistic operating conditions. However, one can also realize an informed projection (i.e., based on local weather forecast) of T inlet,dis in Equation (14). Figure 6 shows simulation results for different set-point temperatures and control horizon setup.  Figure 7 illustrates two additional cases in which some valves are permanently closed with a target outlet temperature of 600 • C and a similar controller setup previously described. The first case (Figure 7a,b) corresponds to valves 2, 3 and 4 being closed, while the second corresponds to valves 1, 2, 4 being closed (Figure 7c,d). The mass flow rates have been plotted during a complete charge-discharge cycle once the steady state is achieved. The objective of these experiments is to show either the effect of inoperable valves due to technical reasons (maintenance or failure), or that of a more affordable plant design (less valves and ducts).
A more flexible operation is achieved in the first case by means of an available valve close to the TES inlet (during discharge operation), for this region of colder rock temperature reduces the outlet temperature more effectively when mixed with that from the TES outlet (valve no. 5). Figure 7b shows the mass flow rate for each extraction during discharge. Due to the inverse numeration of the extraction flows during the charge, the flow of the first discharge valve (ṁ 1 ) equals the total flow during the charge phase (62 kg/s). An almost continuous mix of flows allows the system to operate close to the objective temperature. In the second case (Figure 7c,d), the proposed valve availability limits the ability of the controller to track the desired temperature. It must be acknowledged that the absence of optimal control of the charge leads to this particular outcome.

Conclusions
This work has introduced a novel strategy to optimally discharge segmented PBTES by means of intermediate extractions. Flow multi-extraction provides a means to actively control the outlet temperature during discharge, avoiding the typical temperature degradation in unidirectional systems. Using an NLMPC scheme to track a desired storage outlet temperature, the extraction mass flowrates can be also optimized. This information is then used to solve for the position of the valves of a more detailed pressure-flow model of the plant. Several scenarios were simulated to show the added flexibility of a particular storage design-each scenario presents a different control setup aiming at a different scheduled discharge temperature. Two additional scenarios were designed to evaluate the importance of the position and number of extractions-the expected outcome shows that, as fewer valves are available, those located on the cooler side of the storage become more important to stabilize the temperature of the mix. The successful implementation of NLMPC for tracking showed the potential benefit of multi-extraction in segmented packed beds for power generation, improving the flexibility and reducing the size of the storage at the expense of a more complex airflow system.
The question of the height-to-diameter limit beyond which the radial flows (2D flow) hinder the applicability of the present formulation in non-detached segmented PBTES is left open. For smaller cross-sections the effect of the extractions on the radial temperature gradient may be small, but for larger sections, the assumption of small radial gradients may only be valid if the extraction flow is drawn between physically detached segments of the TES. In the case of large radial temperature gradients, the optimal control approach is still valid as long as a suitable, accurate predictive model is found that can be solved in a reasonable time.
Author Contributions: A.R. developed the control strategy, elaborated the numerical models, and wrote the manuscript. R.C. contributed in the modeling phase, and reviewed the manuscript. E.G. was responsible for securing the funds and reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

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

Abbreviations
Greek Symbols