ORC Optimal Design through Clusterization for Waste Heat Recovery in Anaerobic Digestion Plants

: Waste heat recovery (WHR) systems through organic rankine cycles (ORCs) in anaerobic digestion plants may improve cogeneration efﬁciency. Cogeneration unit power output, ﬂue gas temperature, and mass ﬂow rate are not constant during the day, and the thermal load requested by digesters shows seasonal variations. For this reason, a proper design of the ORC is required. In this study, a design methodology is proposed, based on the clustering of the boundary conditions expected during one year of operation and the anaerobic digestion plant operation. The design has to be a compromise between part-load operation and nominal power rating. In this study, the ORC design boundary conditions were partitioned into four representative clusters with a different population, and the centroid of each cluster was assumed as a potential representative boundary condition for the cycle design. Four different ORC designs, one for each cluster, were deﬁned through an optimization problem that maximized the cycle net power output. ORC designs were compared to those resulting from the seasonal average boundary conditions. The comparison was made based on the ORC off-design performance. Part-load behavior was estimated by implementing a sliding-pressure control strategy and the annual production was therefore calculated. ORC off-design was studied through a detailed Aspen HYSYS simulation. Simulations showed that the power output of each design was directly connected to the cluster population. The design obtained from the most populated cluster generated 10% more energy than that from a system designed by taking into account only the year average conditions.


Introduction
Anaerobic digestion plants have been spreading in Europe due to the favorable synergy between energy production and sanitation services [1]. Moreover, subsidies related to biogas upgrading into biomethane and its injection in the natural gas network, or its liquefaction, have made anaerobic plants profitable in several European countries [2]. All these issues make anaerobic digestion a promising waste-processing technology from environmental, energy, and economic points of view [3]. Anaerobic digestion requires electric energy for sludge piping and handling and thermal energy for maintaining the correct digestion temperature [4]. Traditionally, biogas-fueled internal combustion engines (ICEs) are used as cogeneration units [5]. However, micro-gas turbines (mGT) are preferred in some cases due to the lower maintenance cost and pollutant formation [6]. The thermal energy produced by cogeneration units is often more than that required by digesters, especially in hot climates or when a mGT is adopted. The excess thermal energy may be exploited in a waste heat recovery (WHR) system [7]. Several WHR systems have been proposed to increase the efficiency of biogas plants. Bruno et al. [8] proposed mGT inlet air cooling through an absorption chiller that exploits the excess thermal energy in the flue gas. Results showed that for anaerobic digestion plants equipped with a mGT, inlet air cooling is a valid technology from energy and economic points of view. A further study [9] showed that inlet air cooling effectiveness is strongly dependent on the climate conditions. In [7,10], it was demonstrated that biogas plants could be optimized to become a polygenerative system from which cooling, freshwater, and electrical power are obtained. In a previous study by the authors [11], the biogas plant capability of powering both district heating and district cooling was analyzed. Results showed that trigeneration was a viable improvement for an anaerobic digestion system, but the small thermal load required in mid-seasons penalized the payback period.
Among all the possible WHR solutions for power production, organic rankine cycle (ORC) is the most diffused and studied. The simplicity [12], flexibility [13], and technological maturity [14] of this technology have led to the proliferation of ORC modules in several countries [15].
In the last few years, several ORCs have been installed in anaerobic digestion plants to increase electricity production. Many of them have come from downsizing and customization of commercial modules [16,17]. In a previous study [18], the authors showed that ORC is a promising technology to increase anaerobic digestion plant efficiency, even by adopting commercial modules. In [19], the authors showed that even small modifications of an anaerobic digestion plant could increase the producibility of a commercial ORC and the global efficiency of the cogeneration system.
Many authors have focused their research on the optimization of ORC modules for biogas applications. In [20], Dumont et al. focused on analyzing three different cycle architectures: subcritical, wet-expansion, and supercritical ORC. The analysis considered fixed design conditions for the anaerobic digestion plant and the results showed that the supercritical cycle maximized the investment net present value (NPV). In [21], Mudasar et al. proposed an ORC coupled to a biogas combustor and chose the best configuration for fifteen different operating conditions, which might occur due to changes in biogas composition. The work did not consider the ORC module part-load or off-design operations. In [22], Benato et al. determined the optimum plant configuration and the optimum ORC operating fluid among 118 candidates to maximize an ORC performance coupled to a biogas engine. In [23], Sevinchan et al. presented a system for WHR from a biogas fueled mGT by using both an ORC and an absorption chiller. The system achieved an energy efficiency of 72.5% and a second law efficiency of 30.44%. In [24], Hosseini et al. proposed the optimization at nominal conditions for an ORC powered by the exhaust gases of a biogas fueled mGT. Several designs were proposed for different mGT operating conditions. In [25], Linnemann et al. studied an ORC for biogas engine WHR by optimizing the evaporator geometry. All the studies above-mentioned considered fixed conditions as boundary conditions for the ORC design, and they neglected any variations in the anaerobic digestion plant operation. However, anaerobic digesters may change their thermal requirements daily, with significant seasonal variations [26], thus profoundly influencing ORC performance.
Moreover, ambient temperature variation during the day influences the exhaust gas temperature of a mGT cogeneration unit and therefore the ORC evaporation temperature [11]. These phenomena should be considered during the ORC design to assess its actual performance. In [27], Sung et al. partly considered this issue by proposing different sizes for an ORC + mGT system by considering the actual plant off-design in twelve different average conditions, one for each month of the year. The system was optimized, and the off-design of both mGT and ORC was analyzed to determine the best configuration. The analysis considered the seasonal modification of the digester thermal load, but did not consider the daily variations of cogenerator performance.
This study proposed a methodology to fill the literature gap regarding the design of an ORC as a WHR unit for a biogas-fueled cogeneration unit. The design was performed by considering the digester load, the varying ambient conditions, and cogeneration unit operation. As a case study, the operating data presented in [11] for an anaerobic diges-tion plant near Pisa (central Italy) were considered. The ORC design was performed by clustering the plant data in four representative clusters and by defining four different designs, one for each cluster centroid. The ORC designs were compared among each other and with the designs related to the year average, winter average, and summer average operating conditions. The objective of the comparison was to select the design with the highest annual production. The comparison purpose was to demonstrate that careful choice of boundary conditions leads to an improved ORC design by maximizing the energy produced in the year. In particular, the energy production increase in comparison to the choice of annual average operating conditions was assessed to quantify the advantage of using the proposed design strategy based on a clusterization procedure.
In theory, ORC design and off-design performance may be jointly considered, and they could be optimized together to maximize the annual production. In practice, this may result in a complex problem from an optimization point of view. The optimization would directly search for a trade-off between ORC power rating and operating hours since both parameters should be maximized. The reason for this trade-off is that a too-large power rating may yield a low production due to part-load operation. However, the search for the optimal configuration from both design and off-design operating conditions may be approximated by carefully choosing the design boundary conditions. The design can be divided into two steps: • a set of representative boundary conditions is selected for the design optimization problem through clusterization.

•
The design optimization problem is solved with a standard approach: the ORC power output is maximized for the given set of boundary conditions. As a previous study demonstrated, this may be efficiently and reliably done with general-purpose solvers [28].
After this, the ORC off-design performance was simulated, and the yearly production was calculated to verify the proposed methodology. A detailed off-design model was created to do this by using Aspen HYSYS, a state-of-the-art software for off-design system simulations.
In summary, the paper aims at demonstrating that data clustering is a useful strategy to determine the design boundary conditions for maximizing the annual energy production, but in a simplified, yet effective way.

Reference Layout
This work refers to an existing anaerobic digestion plant located near Pisa (central Italy), whose configuration is presented in Figure 1.
The plant is based on two anaerobic digesters with a total capacity of 4600 m 3 . The plant is operated for the co-digestion of municipal sewage from the wastewater treatment plant and the organic waste from the municipality collection (referred to as sludge). The nominal plant capacity is 10.8 t/h of sludge. Digesters operate at a constant temperature of 37 • C (mesophilic digestion), producing biogas with a final methane concentration of 65%. The biogas is desulfated, stored in a gasometer, and then used to power a Capstone CR600 microturbine [29], designed to work with biogas, consisting of three modules of 200 KW el each, for a total of 600 KW el .
The produced electrical power meets the internal plant need. The turbine exhaust gases are sent to a heat exchanger, where they transfer heat to an intermediate water circuit that heats the incoming sewage. The thermal energy in excess is dissipated by diverting part of the mass flow before the heat exchanger. Before entering the heat exchanger, a part of sludge inside the reactors is recirculated and mixed with incoming sewage to avoid overheating. The recirculation ratio, defined as the ratio between the mass flow rate of sludge recirculated and the mass flow rate of sludge entering the plant, was around 23:1. The high recirculation ratio allowed for the sludge temperature in the heating system to be kept near 37 • C without damaging the digestion substrate. Reference anaerobic digestion biogas facility layout. For further details, the reader should refer to [11].
The plant behavior was simulated in a previous study through a hybrid steady-state and dynamic model [11]. From the model, operating data on a fifteen-minute basis for a period of one year are available, comprehending, among the others: • Exhaust gases temperature. • Exhaust gases composition, assumed to be constant over time (0.038% CO 2 , 0.020% H 2 O, 0.190% O 2 , and 0.752% N 2 ). • Thermal power required by incoming sewage to keep the optimal temperature in the reactors.
Historical series of climate data for the definition of ambient conditions were obtained from [30]. Figure 2 reports the histograms of plant operating conditions. As it can be noted, there are some preferential operating conditions, especially concerning the heat flow rate required by the digester . Q dig and the mass flow rate of turbine exhaust gas . m f . On the other hand, the turbine exhaust temperature T f resulted in having a more limited variation between its maximum and minimum values.
m f , and T f strongly impact ORC operation. As these quantities significantly vary during the year, a design approach that allows the mitigation of ORC part-load is further justified.

Modified Layout
The study aims at analyzing a new system configuration to improve plant performance and profit. The proposed solution implies a further exploitation of the thermal content of the microturbine exhaust gases. The exhaust gasses are used to heat the sewage, as in the original configuration, and power the ORC. Figure 3 shows the modified system layout.
mGT exhaust gases provide the heat to the fluid of an intermediate circuit, which, conversely to the original scenario, powers the ORC evaporator and then heats the incoming sewage. The incoming sewage is preheated by the sewage exiting the digesters, which is mostly recirculated into the reactors. This internal regeneration reduces the heat requirement of the digesters, even though only to a small extent, and allows a slight increase in the available thermal energy to the ORC. The sewage heating is regulated through an air cooler inserted into the intermediate circuit, upstream of the circulation pump. This equipment rejects the excess of heat to the environment. Compared to a bypass, this solution may allow the heat exchanger of sewage heating to work at a constant flow rate under optimal heat-transfer conditions.  Diathermic oil is used in the intermediate circuit to match the operating temperatures. In the analyzed case, Therminol 66, a widely used synthetic fluid, was selected. Such a thermal oil can operate in a vast temperature range (from around −3 to 345 • C) with high thermal stability and long operating life [31].
Concerning the ORC, a subcritical, internally regenerated configuration was chosen. Internal regeneration might be convenient in the analyzed case, as the mGT waste heat is used for the ORC and meets the digester thermal demand.
Pentafluoropropane (R245fa) was selected as the ORC working fluid. In the literature, R245fa is referred to as one of the most appropriate fluids for medium-low temperature ORC applications, from both economic and performance points of view [32]. Furthermore, R245fa has quite good safety characteristic and optimal flammability limits. The main characteristics of the selected fluid are shown in Table 1 ORC cooling is obtained through a water-cooled condenser using purified wastewater from the purification plant. In this way, the auxiliary consumption associated with air coolers is saved, thus allowing for more stable condensation conditions. The optimal design of an ORC may vary from case to case according to the thermal source and other external conditions. In this study, the preliminary plant design was addressed as a constrained optimization problem. A mathematical model of the system was defined and the resulting optimization problem was solved by maximizing the net electrical power produced by the ORC.
The design of the system must be performed according to some boundary conditions, which represent the operational conditions of the mGT and digesters. The ORC design must consider the mGT exhaust gases' mass flow rate . m f and maximum temperature T f and the heat flow required by digesters . Q dig to keep constant the digestion process temperature. These boundary conditions constrain the operation of the ORC since its introduction should modify neither the mGT operating conditions nor those of the digester. On one hand, ORC should be designed according to anaerobic plant operation; on the other hand, the design should aim at minimizing the ORC off-design to maximize the energy recovery over the year. For this reason, the ORC design boundary conditions have to be carefully chosen and should correspond to the working conditions that are more likely to occur during the year. To identify these most representative operating conditions, . m f , T f , . Q dig were clustered. This clustering allows for the identification of the plant's most likely working conditions by considering the correlation between these three variables.
Clustering consists of defining distinct groups (clusters) that are considered representative of data and assigning data to these groups. The distribution is carried out in such a way that the data exhibited substantial similarity within the same group and strong dissimilarity with the elements of other clusters. In the current study, the k-means algorithm was used, as implemented in MATLAB v. 2019b [35], according to Lloyd's algorithm [36]. The k-means algorithm is the most representative partitioning algorithm. In these algorithms, cluster membership is defined based on the distance of the element from a representative cluster point (centroid). Other clustering techniques are available, but if a considerable number of data are available, partitioning techniques are more suitable as they are more straightforward and efficient.
As the k-means algorithm defines centroids to minimize the sum of the distances with data, the underlying problem is a non-linear optimization problem. As it often happens, these problems tend to be trapped into local minima. A multi-start approach has been used to avoid this, and the clustering algorithm was started from several different starting points, and the best solution at the end of all the runs was chosen. This technique made the analysis repeatable and helped to avoid suboptimal classification. The least number of random starts that led to stable results was found equal to 10 through a trial and error approach.
For the k-means algorithm, the number of centroids must be specified in advance. This requirement might lead to non-natural partitioning due to the non-representative number of partitions. When a low number of clusters is specified, two or more different groups might gather together, yielding a misleading partitioning. The same can happen if there are too many clusters: elements belonging to the same group might be erroneously partitioned into different groups. Partitions from two to 10 clusters were tested to determine the optimal number of clusters for the available dataset. The average silhouette value was used Appl. Sci. 2021, 11, 2762 7 of 24 to measure the partitioning goodness. The silhouette is a standard indicator to understand if the right number of clusters is used and is defined as in Equation (1) [37]: where a i is the average distance between the i-th point and the other points in the same cluster as i, while b i is the minimum average distance between the i-th point and the points of a different cluster, minimized between the different clusters. Silhouette values range between −1 and 1, and for values close to 1, it indicates the correct assignment of the point to the corresponding cluster. A null silhouette value indicates that the point does not belong in a defined manner to one cluster compared to another, while a negative value indicates that the point is assigned to the wrong cluster. Clustering results are reported in Figure 4a, which show how the four-cluster case performed the best according to average silhouette values. The silhouette value distribution across the four-clusters case is shown in Figure 4b. Since several data points within the same clusters scored high silhouette values, it can be concluded that the partitioning was successful.
starting points, and the best solution at the end of all the runs was chosen. This techniqu made the analysis repeatable and helped to avoid suboptimal classification. The leas number of random starts that led to stable results was found equal to 10 through a tria and error approach. For the k-means algorithm, the number of centroids must be specified in advance This requirement might lead to non-natural partitioning due to the non-representativ number of partitions. When a low number of clusters is specified, two or more differen groups might gather together, yielding a misleading partitioning. The same can happen there are too many clusters: elements belonging to the same group might be erroneousl partitioned into different groups. Partitions from two to 10 clusters were tested to deter mine the optimal number of clusters for the available dataset. The average silhouette valu was used to measure the partitioning goodness. The silhouette is a standard indicator t understand if the right number of clusters is used and is defined as in Equation (1)  (1 where ai is the average distance between the i-th point and the other points in the sam cluster as i, while bi is the minimum average distance between the i-th point and the point of a different cluster, minimized between the different clusters. Silhouette values rang between −1 and 1, and for values close to 1, it indicates the correct assignment of the poin to the corresponding cluster. A null silhouette value indicates that the point does not be long in a defined manner to one cluster compared to another, while a negative value in dicates that the point is assigned to the wrong cluster. Clustering results are reported i Figure 4a, which show how the four-cluster case performed the best according to averag silhouette values. The silhouette value distribution across the four-clusters case is show in Figure 4b. Since several data points within the same clusters scored high silhouett values, it can be concluded that the partitioning was successful.  Due to the inherent noise in the large amount of data available, clustering was performed on hourly averaged data, rather than on the entire dataset. This timescale reduction allowed for the achievement of more stable results. This design choice is considered as adequate if the slow dynamic of the anaerobic digester is taken into account. Table 2 shows the cluster centroid values obtained from clustering. The values chosen to size the system are those corresponding to the most populated cluster (Cluster 3). These values also correspond to those capable of maximizing the ORC cycle electrical power output. In Table 2, the four clusters are compared with the annual average plant operating conditions, and with the summer and winter average conditions. However, it may be noted that the clustering algorithm allowed for the selection of representative operating conditions, which cannot be selected with traditional methods. Data clustering may disclose a pattern in data that would be difficult to highlight otherwise. For example, in Table 2, the system operating conditions that corresponded to the four cluster centroids were compared with the annual and summer average conditions. As a result, the annual and summer average did not resemble any of the clusters, which means that the former failed in representing the actual frequency of the operating conditions. In other words, the annual or summer average was not like the operating conditions that frequently showed up in the record, while the cluster centroids were precisely built to resemble those that most frequently showed up. Therefore, each cluster represents a part of the operating condition dataset and the number of elements represented by each centroid (i.e., similar to it) is the cluster population.
As a result, Cluster 3 represents most of the operating conditions over the year. This means that this cluster is similar to the operating conditions encountered most frequently during system operation.
Clusters are built through a numeric procedure, so they are mathematical objects. However, they may maintain physical significance. In this case, the following patterns can be observed: To summarize what has been done thus far, some boundary conditions that represent the operational condition of the plant are needed to perform the ORC design. A clustering technique (k-means) was used for defining a representative set of these conditions. A four-cluster partition described the whole dataset in the best way. Among the four clusters, the most populated one, which accounts roughly for a third of the operational conditions, was selected as the boundary condition for the ORC design. Based on such values of . m f , T f , and . Q dig , the ORC was optimized to maximize its electrical power output.

Design Optimization Problem
The preliminary ORC design was formulated as a constrained optimization problem, which was solved through MATLAB v. 2019b.
CoolProp v. 6.3.0 [38] was used to determine the thermodynamic properties of the ORC fluid. In the mathematical model, each component is assumed to be an open system operating in steady-state conditions, and the pressure drops in ORC heat exchangers and pipes are neglected.
The part of the digestion plant modeled in the code includes the ORC and the thermal fluid loop ( Figure 5). The design is performed by considering that all the available thermal energy is exploited under design conditions. Therefore, in this preliminary stage, the presence of the air cooler located upstream of the thermal oil circulation pump is overlooked. This heat exchanger is used to reject the excess heat to prevent sewage overheating, and its operation is considered only in off-design conditions. The ORC design problem can be formulated as in Equation (2): where Φ is the so-called feasible region; f (x) is the problem objective function; and x is the variable vector. The optimization problem must satisfy the following set of constraints (Equation (3)): where lb and ub represent the lower and upper variable bounds, respectively; A and b define the linear inequality set that must be observed by x components; finally, c(x) is a vector function so that c: R 6 → R n , where n is the number of the problem's non-linear constraints. c(x) ≤ 0 represents a set of non-linear inequalities that the problem must satisfy. The x components (i.e., optimization variables) are: Thermal oil temperature exiting the ORC evaporator T loop,2 in • C. For thermal oil loop variable subscripts, please refer to Figure 5. The ORC cycle was modeled to evaluate f (x) in Equation (2). The following equations define the calculation procedure for each ORC point. This equation set must be interpreted by considering that the optimization algorithm tries a set of optimization variables (x components) in each solver iteration. Given that during each ORC cycle evaluation, T ev,orc , T cd,orc , and ∆T sh,orc are known, among the other variables, since the solver provides them. Since the ORC evaporation temperature is known, it is possible to determine ORC evaporation pressure p ev,orc by using CoolProp. The ORC turbine inlet conditions are then univocally determined (p = p ev,orc ; T = T ev,orc + ∆T sh,orc ). Since T cd,orc is known, condensation pressure p cd,orc may be readily determined. At this point, isentropic expansion condition may be calculated and, by following the definition of turbine isentropic efficiency η is,exp , (equal to 0.8), the turbine specific work ∆h exp is determined (Equation (4)): For the interpretation of ORC-related subscripts, please refer to Figure 5.
Depending on the condensation pressure and the temperature at the pump inlet (T cd,orc − ∆T sc,orc ), pump inlet enthalpy and entropy were obtained. By knowing the evaporation pressure, the pump outlet isentropic conditions can be estimated. Through the pump isentropic efficiency η is,pump , equal to 0.7, the pump specific work ∆h pump is calculated (Equation (5)): Working fluid mass flow is determined from the evaporator thermal input . Q ev and evaporator net enthalpy difference ∆h ev . The first one is calculated from T loop,1 , T loop,2 , and . m loop , which are provided by the solver in each iteration (Equation (6)): ∆h ev is calculated by determining the regenerator hot side outlet conditions. These are calculated from the regenerator thermal effectiveness ε reg defined in Equation (7): Regenerator hot outlet conditions may be used to define regenerator cold side outlet enthalpy h 2r,orc (i.e., evaporator cold side inlet enthalpy) through regenerator energy balance, as reported in Equation (8): From this, the evaporator net enthalpy difference can be calculated (Equation (9)): Finally, by combining Equation (6) and (9), . m orc is calculated as in Equation (10): From Equations (4), (5), and (10), ORC net power output W net,orc is determined as in Equation (11): where η el is the electro-mechanical conversion efficiency, taken as equal to 0.95. Theoretically speaking, W net,orc might be used as the optimization problem objective function (f (x) in Equation (2)). Although this allows for the power output global maximum determination, it was noted that the same (maximum) W net,orc value is yielded by several different cycle configurations. To univocally identify the best solution, the objective function was modified by including the thermal oil loop circulation pump electric consumption, defined as in Equation (12). The computation of the pump consumption in the loop allows for the optimum to be reached by a unique set of variables. where: • η is,pump,loop is the oil circulation pump isentropic efficiency (η is,pump,loop = 0.7). • ∆p loss,i is the estimation of the pressure losses in each i-th exchanger of the loop (10 kPa). • ρ loop is the thermal oil density at the pump inlet.

•
The requirement of considering the pumping power is more related to the solver numerical behavior than to the system actual physics. The optimization algorithm may consider as equivalent two solutions that are only very similar. To help it to differentiate between them and since the threshold used to distinguish between the solutions is a very small number, it is sufficient to include the pumping power in the calculations to consistently achieve the same optimum every time, for the same set of boundary conditions. However, despite the significance for the optimization algorithm, the selection of the pressure drop value can be quite arbitrary due to its impact on performance being almost negligible.
The modified objective function is equal to (Equation (13)): The feasible region Φ introduced in Equations (2) and (3) is defined through the optimization problem constraints. For each of the optimization variables, lower bounds (lb) and upper bounds (ub) define the search domain as Equations (14)- (19): T cool,out ≤ T cd,orc ≤ T source,in − ∆T pp,source where T crit,orc is the ORC fluid critical temperature; T cool,out is the cooling water temperature exiting the ORC condenser, which is set to 25 • C; T source,in is the temperature of the flue gases at the inlet of the waste heat recovery heat exchanger; and ∆T pp,source is the minimum allowable pinch point between the flue gases and the thermal oil. The ORC operational limitations were imposed through linear and non-linear constraints. A minimum allowable pinch point was set for each heat exchanger. Both ∆T pp,source (pinch point between exhaust gases and thermal oil) and ∆T pp,sew (pinch point between thermal oil and sewage) were equal to 30 • C. ∆T pp,ev (ORC evaporator pinch point) was equal to 15 • C. ∆T pp,reg (ORC regenerator pinch point) was equal to 10 • C. Finally, ∆T pp,cd (ORC condenser pinch point) was equal to 10 • C.
Maximum ORC temperature T max,orc (turbine inlet) must be at least 5 • C lower than that of the organic fluid degradation. The ORC condensation pressure p cd,orc must be higher than 0.5 bars. A minimum superheating degree of 5 • C was set to prevent liquid injection into the ORC turbine. Similarly, to avoid pump cavitation, a subcooling ∆T sc,orc of 5 • C was set. Other considerations are related to the thermal oil loop, for which a minimum design temperature T loop,min of 70 • C was imposed due to the high oil viscosity at low temperatures. For exhaust gases, a maximum value for the oil heat exchanger effectiveness ε source of 0.8 was imposed to reflect the performance of commercially available heat exchangers. Finally, a minimum mGT exhaust gas temperature T source,out of 70 • C was set to prevent acid condensation. All the previous considerations (and some others) can be translated into constraints, both linear and non-linear. As for the linear constraints, they can be defined as in Equation (20)- (26): T ev,orc + ∆T sh,orc ≤ T loop,1 − ∆T pp,ev T ev,orc + ∆T sh,orc ≤ T max,orc (22) T cd,orc ≥ T cool,in + ∆T sc,orc + ∆T pp,cd (23) T loop,2 ≥ T sew,out + ∆T pp,sew where: • ∆T loop,min is the minimum achievable temperature drop of thermal oil across the ORC evaporator. • T sew,out is the sewage temperature at the outlet of the sewage heating heat exchanger, set at 37.2 • C to consider the thermal losses before entering the digester.
As for the non-linear constraints, the equation c(x) ≤ 0 can be translated into a nonlinear inequality set as defined in Equation (27)-(38): T 2r,orc + ∆T pp,reg ≤ T 4,orc (31) T cool + ∆T pp,cd ≤ T cd,orc T source,out ≥ 70 (34) T sew,in + ∆T pp,sew ≤ T loop,3 where: • χ exp,out is the vapor quality at the end of the expansion. • ∆T reg,min is the minimum temperature difference between the temperature at the end of the regeneration T 4r,orc and ORC condensation temperature T cd,orc . • T 2,orc and T 2r,orc are the inlet and outlet temperatures on the regenerator cold side (evaporation pressure), respectively, while T 4,orc and T 4r,orc are the inlet and outlet temperatures on the regenerator hot side (condensation pressure), respectively. • T loop is the oil temperature at the beginning of organic fluid evaporation (thermal oil temperature at evaporator pinch point). • T cool is the cooling water temperature at the beginning of the organic fluid condensation (thermal oil temperature at the condenser pinch point). • T loop,5 is the thermal oil temperature at the circulation pump outlet. • T loop,3 is the thermal oil temperature entering the air cooler. Under design conditions it is equal to T loop,4 (i.e., the temperature of the oil exiting from the air cooler). • T sew,in is the sewage temperature at the inlet of the sewage heating heat exchanger set at 37 • C, which is the digestion process's nominal temperature.
Some of the quantities involved in non-linear constraint evaluation must be calculated, as they are not readily available. T loop , which is necessary to set the constraint on the evaporator pinch point, was found as in Equation (39): where h 2 ,orc is the organic fluid enthalpy at the beginning evaporation. T cool , which is necessary for the constraint on the condenser pinch point, was derived from Equation (40): h 5,orc is the organic fluid enthalpy at the beginning condensation and T cool,in is the cooling water temperature entering the ORC condenser, which was set to 15 • C. Finally, . Q cond is the condenser heat flux, which was found as (Equation (41) m f and T f ), previously calculated through clustering, T source,out was calculated as in Equation (42): Likewise, T loop,5 , which is necessary to constrain the pinch point of the flue gases in the heat exchanger, was determined as in Equation (43): Thermal oil temperature exiting the sewage heat exchanger was calculated to impose the oil loop minimum temperature: Finally, T sew,in was calculated to set the constraint on the pinch point of the sewage heat exchanger: where . m sew and c p,sew correspond to sewage mass flow rate and specific constant pressure thermal capacity, respectively.
To solve the problem defined by Equations (2) and (3), a sequential quadratic programming (SQP) algorithm was used, as implemented in [39], inspired by the SQP algorithm presented in [40]. The problem is not smooth in theory, as the use of CoolProp to evaluate fluid thermophysical properties makes the objective function a "black box". In this case, the objective function may be numerically estimated, but it does not have derivatives. Despite this, the problem is sufficiently smooth to be solved with an algorithm based on derivatives. The numerical estimation of gradients and Hessians is most often successful due to the general smoothness of thermophysical properties.
The resulting optimization problem is non-linear, without an analytical representation, and potentially non-convex. In this case, the solver might find a local minimum instead of a global one. Such an issue may be avoided by repeating the optimization from several starting points. A multi-start algorithm has been applied to improve the convergence rate to optimum, as implemented in [41,42].

Preliminary Design Results
Once that the whole problem has been implemented, it can be solved for the operational case determined by the clustering reported in Section 3.1.1. The optimization variable values that correspond to the most populated cluster (Cluster 3) are reported in Table 3 as well as the objective function (ORC net power output) value.

Off-Design
Once the design conditions were known, an off-design model was developed in Aspen HYSYS. Aspen HYSYS is a simulator used to model chemical reactive and non-reactive processes from unit operations to full plants. For a given design condition, the software can design and simulate the behavior of various typologies of commercial heat exchangers such as shell and tube, plate multi-stream heat, and air coolers.
The off-design model of the ORC was implemented by considering design data from MATLAB optimization. The model considered the main components of the ORC, thus including heat transfer fluid loop and digesters. Shell and tube heat exchangers were considered for the evaporator, regenerator, and water-cooled condenser. Each heat exchanger was designed in Aspen Exchanger Design and Rating (Aspen EDR) and then imported in Aspen HYSYS. Starting from the geometry and the flow conditions in the heat exchanger, the model created an internal mesh to solve the heat transfer problem (finite volume approach). In each volume, the code estimated the actual value of the overall heat transfer coefficient (U) by adopting Aspen heat transfer and fluid service (HTFS) correlations [43]. By repeating the calculation in each finite volume, the heat transfer process was solved, and the flow conditions at the heat exchanger outlet were determined.
The solver determined the oil pump and the ORC pump off-design behavior by considering the design conditions and the expected head in three operating points. Commercially available data were adopted to determine the operative curves. Centrifugal pumps were supposed for both the ORC and the heat transfer loop. Regarding the turbine, a more complex criterion based on Stodola's law was applied to determine off-design behavior, since the software did not provide any curve for commercial turbines. The approach followed to derive the turbine off-design curves is detailed in the next subsection.
Digesters were simulated as a heat sink due to the lower operating temperature, which did not provide any critical issue in the heat transfer process.
As for ORC management, a sliding-pressure control strategy was assumed by optimizing the evaporation pressure, superheating grade, and mass flow rate for each operative point. Optimization was performed in Aspen HYSYS by maximizing the system net power output.
The result of the ORC module off-design analysis was an operating map, which is function of the mass flow rate, the temperature of the flue gas, and the heat requested by the digesters. From the annual turbine and digester behavior, the actual energy produced by the ORC module was estimated by interpolating the obtained maps.

Turbine Off-Design Analysis
Stodola's ellipse law provides a simplified approach for predicting off-design performance in the idealized case of an infinite stage number turbines, where acoustic choking could never happen [44]. In the case described in this paper, due to the low evaporating temperature obtained, a single-stage turbine was enough to avoid flow choking. Stodola's approach may be modified to account for a finite number of stages, up to single-stage turbines, as demonstrated in [44]. Stodola's approach is based on the idea of assimilating an entire group of stages to a single nozzle, for which the following proportionality relation holds (Equation (46)).
where B is the turbine outlet static pressure; p 0 is the total inlet pressure; and ψ is the flow coefficient, defined as in Equation (47): where . m is the mass flow rate through the turbine and v the specific volume at the turbine section at the pressure p 0 . Equation (46) is called Stodola's ellipse as it may be represented as an ellipse in the plane.
In [44], it was demonstrated how the ratio ψ * (0 ≤ ψ * ≤ 1) between ψ in any condition and ψ at acoustic chocking condition may be written according to Equation (48) for a finite number of stages: where a is the acoustic chocking expansion ratio. To maintain physical significance, ψ * follows Equation (48) as long as B > p 0 ·a. For B = p 0 ·a, ψ * reaches its maximum value (i.e., ψ * = 1) and for lower B values (i.e., beyond acoustic choking), ψ * = 1 is maintained as the flow coefficient does not vary. By following Equation (48), the proportionality relation in Equation (46) may be restated as in Equation (49) [44]: where the d subscript identifies design quantities. From Equation (49), the off-design relations for finite stage number turbines may be derived as in Equations (50) and (51) for B > p 0 ·a and B ≤ p 0 ·a, respectively [44]: Appl. Sci. 2021, 11, 2762 16 of 24 and: Equation (29), which is valid for chocked operations, stems from the fact that ψ * and ψ are constant for B ≤ p 0 ·a. Both Equation (50) and (51) link off-design pressures to design specifications and allows for the determination of off-design operating conditions.
In the case of real gases, like ORC working fluids, a must be calculated numerically, as the ideal gas approach does not provide a good estimation of this parameter. Concerning the ORC design provided in the previous section, the expansion ratio is well below the acoustic chocking one, even considering strong off-design conditions. Therefore, in the following analysis, Equation (51) will be used, instead of Equation (50), which is used for non-chocked conditions.
In addition to Equation (51), some additional assumptions must be made to define offdesign operation completely. In this case, turbine outlet pressure and inlet temperature are assumed to be constant and equal to their design values. Of course, different assumptions could be made for different control strategies.
By assuming some inlet and outlet conditions and by exploiting Equation (51), offdesign isentropic enthalpy difference across the turbine (isentropic specific work) may be calculated. However, some further corrections must be considered. As suggested in [45], turbine isentropic efficiency should be correct in off-design conditions by considering two parameters, CF 1 and CF 2 , to modify η is,d . In Equation (52) [45]: where f 1 and f 2 are the functional representation of CF 1 and CF 2 (see [45] for the details).

Off-Design Analysis Results
During the year, the thermal energy available for the ORC varies due to turbine exhaust gas flow rate and outlet temperature fluctuations, and due to variations of the heat flow required by the digesters. These parameters are the ORC design boundary conditions. Since both the ORC and the gas turbine are very quick to adjust their output according to ambient temperature and the digesters' heat requirement, the dynamic system behavior is not a limiting factor. Therefore, the analysis was carried out in the steady-state conditions, and hour time intervals were considered for estimating off-design system operation.
According to sliding pressure control, the ORC adapts to load changes by adjusting the evaporation pressure. The resulting superheating grade and working fluid mass flow rate were determined to maximize the cycle power output. In each operating condition, the result is an optimized cycle, like those reported in Figure 6. The ORC is considered turned off when the net power output is less than 20% of its nominal value.
A grid of operating conditions was mapped and interpolated to define a complete performance operating map such as those reported in Figure 7.

Annual Production
The ORC hourly operation for the year may be estimated from the historical hourly trend of the exhaust gases' mass flow rate and outlet temperature, and by considering the digesters' thermal energy requirement for each hour. Data used for the analysis were the same shown in Section 2.1. ORC net power output histogram is presented in Figure 8 by considering the design derived from the most populated cluster (Cluster 3). ORC operates most of the time in off-design conditions. For around 700 h (less than 1/10 of the year), it operates at nominal power. Furthermore, for over 200 h/y, the ORC cannot operate due to the load being too low (P ≤ 0.2 P rom ). However, the ORC positively affects the plant operation, as it substantially increases the total power output in comparison to the reference configuration. The net power output from the plant is reported in Figure 9. The whole plant net power output is roughly shifted upward by a value equal to the ORC average power output. This value may be calculated by considering ORC off-design operation ( Table 4).    Table 4 shows that the fuel utilization efficiency of the whole plant η tot (mGT + ORC) was increased by a ∆η = 10%, in comparison to the case with no recovery. Furthermore, ORC had a quite high value of equivalent operating hours h eq,orc , which implies a reasonably good use of the installed ORC capacity, despite the off-design and downtime periods. The difference between the mGT exhaust gases' energy content and the biomass digester heat requirement is the gross thermal energy available for the ORC. By dividing this amount for the net ORC production, the annual ORC efficiency η year,orc may be calculated. For the considered design, the average yearly efficiency of the ORC η year,orc was equal to 8%, as reported in Table 4.

Comparison with Alternative Designs
Despite the clear advantages in terms of additional power, the ORC operated in offdesign most of the time. For this reason, the installed capacity might be underused, and an alternative design could be considered. The choice of ORC nominal power output is subjected to a trade-off between installed capacity, which should be maximized to increase the production, and the actual availability of excess thermal energy, given the mGT exhaust gas production and the digester heat requirement. The risk is to overestimate the ORC capacity and to work at low loads. Lower ORC power ratings might characterize alternative designs, which may produce less at full load, but they might be able to operate longer at higher relative loads (i.e., with better efficiencies). This possibility must be ruled out to confirm whether basing the design on the clustering of the ORC operating boundary conditions is a good strategy, or not.
For the comparison with the proposed ORC design, other optimal designs related to alternative sets of boundary conditions were defined. Sensible alternative sets of boundary conditions were represented by the other clusters defined in Section 3.1.1. by the simple annual average, and by two representative seasonal averages (winter and summer). The use of the annual average conditions is a followed approach in literature. Therefore, it is of particular interest to compare the performance of the ORC designed based on the most populated cluster with the one designed based on the annual average operating conditions.
The optimal ORC configurations resulting from the boundary conditions in Table 2 were derived by solving the same optimization problem of Equations (2) and (3). The optimization problem constraints were the same for each case (Table 2), as the solved optimization problem is always the same. The only required adjustment is to update constraints with the new values for the boundary conditions listed in Table 2.
As far as the off-design of such new ORC systems is concerned, the same off-design map of Figure 7 was used. Since the off-design map may be reformulated in terms of variations with respect to design conditions, this map can be used for alternative designs. These designs are centered over different values of . m f , T f , and . Q dig . From a physical point of view, this is a simplifying hypothesis, as the most rigorous approach would require the creation of a new off-design map for each design. However, it might be assumed that the performance of similar systems changes in the same way when their boundary conditions vary. This assumption is supported by the fact that all considered systems are of comparable sizes, and thus they share the same technologies.
The power output histograms as a function of operating hours for the seven considered ORC systems (the proposed one, plus the six alternatives) are provided in Figure 10.  Table 2. (a) The four data clusters. (b) All the other configurations (annual, summer, and winter average) in Table 2.
In Table 2, it may be observed how the configuration related to Cluster 3 (i.e., the most populated one) had the expected highest nominal power output among all the clusters. Cluster 3 was the cluster with the highest mGT exhaust gas mass flow rate and the lowest digester thermal requirement. As also discussed for the summer average conditions, in this case, the ORC had the highest amount of available thermal energy. Accordingly, the optimizer selected a high nominal power for the Cluster 3 ORC, and this was the reason why such ORCs are operated most often in part load, as seen in Figure  10.
However, it is the cumulative energy production in the whole year that measures the design effectiveness. In these terms, as shown in Table 5, the configuration related to Cluster 3 performed better than the other clusters, despite it operating at low loads more often than the other configurations (Figure 10). This result demonstrated that a balance between a high number of operating hours at full load and a continuous off-design operation must be searched. This is especially true for the ORC, which is a flexible component that can operate at several load levels. Clustering actively helps in selecting the best configuration for a trade-off between installed power and operating hours at full load. The ORC design based on the most populated cluster outperformed not only the other clusters, but also the configurations based on the annual average, winter, and summer average conditions. This is a remarkable result since it proves the effectiveness of the design strategy based on clustering. The annual average failed to characterize the variation of plant operating conditions over the year correctly. The design approach based on the annual average, underestimated the ORC production potential, which led to an ORC suboptimal configuration with a too small installed power. The same was true for the winter average, which resulted  Table 2. (a) The four data clusters. (b) All the other configurations (annual, summer, and winter average) in Table 2.
The three configurations related to the discarded clusters were smaller in size than the configuration related to the most populated cluster (Cluster 3). For this reason, they may operate for a larger share of the year at nominal power output. The configuration related to the annual average was similar to that of Cluster 1, whereas the winter average and the summer average yielded the smallest and the largest configurations, respectively. Such a size distribution was due to the different seasonal operating conditions. In winter, the maximum heat flow is required by the digesters, and the minimum thermal energy content is found in the mGT flue gases. Therefore, the heat flow available for the ORC is a minimum, leading to an ORC operating at high equivalent hours, but with a small power output ( Figure 10 and Table 5). In summer, the situation is the opposite. Therefore, a larger ORC nominal power output may be installed. In Table 2, it may be observed how the configuration related to Cluster 3 (i.e., the most populated one) had the expected highest nominal power output among all the clusters. Cluster 3 was the cluster with the highest mGT exhaust gas mass flow rate and the lowest digester thermal requirement. As also discussed for the summer average conditions, in this case, the ORC had the highest amount of available thermal energy. Accordingly, the optimizer selected a high nominal power for the Cluster 3 ORC, and this was the reason why such ORCs are operated most often in part load, as seen in Figure 10.
However, it is the cumulative energy production in the whole year that measures the design effectiveness. In these terms, as shown in Table 5, the configuration related to Cluster 3 performed better than the other clusters, despite it operating at low loads more often than the other configurations ( Figure 10). This result demonstrated that a balance between a high number of operating hours at full load and a continuous offdesign operation must be searched. This is especially true for the ORC, which is a flexible component that can operate at several load levels. Clustering actively helps in selecting the best configuration for a trade-off between installed power and operating hours at full load. The ORC design based on the most populated cluster outperformed not only the other clusters, but also the configurations based on the annual average, winter, and summer average conditions. This is a remarkable result since it proves the effectiveness of the design strategy based on clustering. The annual average failed to characterize the variation of plant operating conditions over the year correctly. The design approach based on the annual average, underestimated the ORC production potential, which led to an ORC suboptimal configuration with a too small installed power. The same was true for the winter average, which resulted in a too conservative design. Winter average conditions led to the smallest ORC installed power, which may be useful only if an ORC design that operates most of the time at full load is looked for. Conversely, summer average conditions yielded the largest installed ORC nominal power. However, from Table 5 and Figure 10, summer average conditions slightly overestimated the ORC potential, and the resulting production was lower than that of Cluster 3. This result proves the ability of clustering in helping to select the optimal ORC configuration. Furthermore, the fact that summer average conditions yield a slightly suboptimal configuration is likely to be dependent on the specific climate zone in which the investigated plant is located. From this standpoint, clustering provides a more robust and generalizable design approach.
In Table 5, the results related to all the investigated configurations can be found. As a design guideline, it can be noted that the ORC energy output followed the population of the clusters: the higher the cluster population, the larger the energy output. In other words, in the case of a great amount of data for ORC design boundary conditions, data clustering may be a useful design tool for WHR applications. This design approach determines the correct design configuration, without requiring the optimization of the plant design and off-design performance at the same time.

Conclusions
In conclusion, the use of an ORC to recover the excess thermal energy of the cogeneration in an anaerobic digestion plant was studied. The part of the plant under investigation (thermal oil loop and ORC unit) was modeled, and its annual operation was simulated. Two models with different detail levels were developed. First, the system was modelled in MATLAB, with a simplified approach based on system constitutive equations and the components' constant performance parameters (isentropic machine efficiencies and heat exchanger effectiveness). Second, the system was modeled in ASPEN HYSYS, which allowed the selection of commercial components (i.e., heat exchangers). Through this detailed model, the part-load performance of the system was simulated. These two models were used to define an ORC design procedure, which considered both the design and the off-design ORC operation.
Theoretically, design and off-design performance could be optimized together to maximize the annual ORC production. However, to simplify the design optimization problem, the design procedure was divided into two steps. The first step was based on the clustering of the ORC boundary conditions (i.e., the biogas facility operating conditions). Through clustering, the most representative (i.e., probable) design boundary conditions may be selected. After this, the ORC design optimization problem can be formulated more straightforwardly (i.e., by maximizing the power output ORC for the fixed set of design boundary conditions).
The biogas facility operating conditions were clustered into four representative clusters. Furthermore, three additional sets of boundary conditions were considered to compare the followed design approach with a more traditional one. Annual, winter, and summer average operating conditions were considered for the comparison. By comparing the performance of the seven alternative ORC designs, it was demonstrated that the design based on the most populated cluster yielded the highest annual production. Such a result demonstrates that using the most populated cluster as design boundary conditions may help to define a design that works well for a larger share of the total cases. Therefore, clustering should be used to select the appropriate design conditions based on which the optimal ORC design was performed. The fact that the use of the most populated cluster may help to maximize the annual production suggests that the clustering should be used as a preliminary step in the design procedure. In this way, the effects due to time-varying operating conditions may be considered in the design process, without explicitly including them in the design optimization problem. This is a simplified approach, but it allows relevant increments in the produced energy. By using the clustering, a produced energy increment of around 10% was achieved in comparison to a design based on the annual average operating conditions.