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

energies Article Modeling, Simulation and Optimal Operation of Multi-Extraction Packed-Bed Thermal Storage Systems † Alberto Romero 1,* , Ricardo Chacartegui 2and Emanuele Garone 1 1Department of Control Engineering and System Analysis, Université Libre de Bruxelles, 165/55 Brussels, Belgium; egar[email protected] 2Departamento de Ingeniería Energética, Universidad de Sevilla, 41004 Sevilla, Spain; ricar[email protected] *Correspondence: ajromerofr[email protected] † This paper is an extended version of our paper published in ECOS 2019. Received: 25 December 2019; Accepted: 17 April 2020; Published: 4 May 2020   Abstract: 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 fluid and rocks as a storage medium. Compared to its alternatives, however, PBTES presents a limited flexibility of operation due to the conventional unidirectional flow, which involves the progressive reduction of the outlet temperature during discharge and thus lowers the thermodynamic efficiency 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 flexibility 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 benefits 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 profiles. 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. Keywords: solar thermal power; segmented packed bed; non-linear model predictive control (NLMPC) 1. 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 – 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 Energies 2020,13, 2247; doi:10.3390/en13092247 www.mdpi.com/journal/energies Energies 2020,13, 2247 2 of 13 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). 2. 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 Energies 2020,13, 2247 3 of 13 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 – 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? 3. Materials and Methods Figure 1shows 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. (a) (b) Figure 1. ( a ) Example of a conical-shaped storage discretization in elements of height ∆x , and (b) discharge of a segmented packed bed storage with five extractions. 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. 3.1. 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 Energies 2020,13, 2247 4 of 13 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]. 3.2. 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 R1 , R2· · · R5 , while the connecting branches are represented as Ra , Rb· · · Re . When the system operates in uniaxial mode, all valves except for R5 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. (a) (b) Figure 2. (a) Network for a 10 elements model with five extractions, and (b) element input/output. ∆pi=pout i−pin i=4.466Rep−0.2ψ0.696ε−2.945e(11.85(log10 ψ)2)(liG2 i)/(ρfdp)(1) Equation (1) describes the relationship between the pressure drop at each element ias a function of the mass flow rate (trough the mass flux G=˙ m/Across ) 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: εi dρihi dt =−1 Vi [(ρuhεA)i−(ρuhεA)i−1]−αv,i(Tf,i−Ts,i)−Ac,i Vi hwall,i(Tf,i−Tinf)(2) h=ZT Tref Cp,fdT (3) (1−εi)ρs des,i dt =−αv,i(Tf,i−Ts,i) + 1 Vikeff dTs dx Ai −keff dTs dx Ai−1(4) Energies 2020,13, 2247 5 of 13 es=ZT Tref CsdT (5) αv,i=700(Gi/dp)0.76 (6) 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. nnode ∑ i=1 ˙ mi=0 (7) Tnode,out −Tref = nnode ∑ i=1 ˙ miCpi(Ti−Tref)/nnode ∑ i=1 ˙ miCpi (8) ∆pbranch,j=RjG2 j(9) Tr(t=0) = Tr,0 (10) T(x=0) = Tin (11) p(x=0, L) = p0,pL(12) ∆˙ m(x=x1,x2, ..., x5) = −[˙ mext1;˙ mext2; ... ˙ mext5](13) Equations (7)–(9) express, respectively, mass–energy and momentum conservation, for every node iand branch jacross 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 ˙ mext,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/m3 , which simplifies the assumption of constant fan pressure so that the prescribed flowrate can be obtained. The specific heat for air was assumed constant ( Cp,f = 1050 J/kg K ), while for the rock a linear interpolation of the data from [ 9 ] was used ( Cp,r = 800 + ( 1100 − 800 )/( 650 − Energies 2020,13, 2247 6 of 13 150 )(Tr− 150 )J/kg K ). Finally, the effective thermal conductivity was approximated by keff = 0.5 + (2.0 −0.5)/(650 −150)(Tr−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]. 3.3. 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 3describes 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: minimize Tf,Tr,˙ mJ= Hc ∑ q=1 γq(Tout,q−Tset)2 subject to Tf,i,k−Tf,i+1,k−ei,k(Tf,i,k−Tr,i,k) = 0 Mr,iCp,r(Tr,i,k+1−Tr,i,k)−ei,k next ∑ q=i ˙ mq,kCp,f(Tf,i,k−Tr,i,k)∆t=0 (Tf,out,k−Tref)Cp,f next ∑ q=1 ˙ mq,k− next ∑ q=1 ˙ mq,kCp,f(Tf,i,k−Tref) = 0 ˙ mmin ≤˙ mi,k≤˙ mmax;∀i=1, . . . , next ∀k=1, . . . , Hc (14) 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 e 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 ( next = 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. Figure 3. Simulation diagram with control scheme. Energies 2020,13, 2247 7 of 13 4. 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. Table 1. Packed bed properties. Variable Symbol Value Unit Particle diameter d0.03 m Porosity ε0.37 - Sphericity ψ0.7 - Storage diameter Dtop 40 m Dbottom 32 m Storage height H25 m Solid density ρr2800 kg/m3 Charge time tc8 h Discharge time td16 h Charge mass flowrate ˙ mc132 kg/s Discharge mass flowrate ˙ md66 kg/s Insulation: wall, bottom, lid Uoverall 0.02 W/m2K Environmental air temperature Tenv 40 ◦C Soil effective temperature Tsoil 20 ◦C Figure 4shows 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. (a) (b) Figure 4. Rock temperature after charge (solid line) and discharge (dashed line) at the end of ( a ) cycle 1 and (b) cycle 30. 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 Energies 2020,13, 2247 8 of 13 of model predictive controllers, which naturally mitigate the associated model inaccuracies by means of the feedback of state measurements [27]. 5. 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 r2 , where r is the volume ratio V/V0 that takes the value 0.1, i.e., the new design is ten times smaller. The new dimensions are, therefore: Dtop = 27.4 m, Dbottom = 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 5shows 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 5also 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. Figure 5. Outlet temperature with new design (r=0.1) compared with original design (r=1). 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 ( Hc=Hp ), 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 Ra...e = [ 0.01, 0.01, 0.01, 0, 0 ] Pa/kg s −1 (Figure 2a). The value of the control resistances R1...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 ( Tinlet,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 Tinlet,dis in Equation (14). Figure 6shows simulation results for different set-point temperatures and control horizon setup. Figure 6a–c represents control horizons ( Hc ) of 12, 14 and 16 h respectively, starting at the end of the charge cycle. In contrast, Figure 6d presents a scenario with Hc = 16 where only for the last 14 h of Energies 2020,13, 2247 9 of 13 the horizon the tracking effort is demanded. Of all four scenarios, only in the first one (Figure 6a) the objective is within the limits of what the controller (its model) and the system can provide. The worst performance is for the scenario with Hc = 16, where the output remains close to a comparatively lower objective temperature (580 ◦ C). Temperature oscillations are significant, the magnitude of which depends on the flow configuration and, more importantly, the control time step (1 h). (a) (b) (c) (d) Figure 6. Output temperature: actual (solid), optimized reference (dash), and uniaxial (dot-dash): ( a ) Tref = 620 ◦ C, 12 h discharge, ( b ) Tref = 620 ◦ C, 14 h discharge, ( c ) Tref = 580 ◦ C, 16 h discharge, and (d)Tref =620 ◦C, last 14 h discharge. Figure 7illustrates 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 ( ˙ m1 ) equals the total flow during the charge phase (62 kg/s). An almost continuous mix of flows allows the system to operate