Analysis of Thermodynamic Models for Simulation and Optimisation of Organic Rankine Cycles

Equations of state (EOSs) form the base of every thermodynamic model used in the design of industrial processes, but little work has been done to evaluate these in the context of such models. This work evaluates 13 EOSs for their accuracy, computational time and robustness when used in an in-house optimisation program that finds the maximum power output of an organic Rankine cycle. The EOSs represent popular choices in the industry, such as the simple cubic EOSs, and more complex EOSs such as the ones based on corresponding state principles (CSP). These results were compared with results from using the Groupe Européen de Recherches Gazières (GERG) EOS, whose error is within experimental uncertainty. It appears that the corresponding state EOSs find a solution to the optimisation problem notably faster than GERG without significant loss of accuracy. A corresponding state method which used the Peng–Robinson EOS to calculate the shape factors and a highly accurate EOS for propane as the reference EOS, was shown to have a total deviation of just 0.6% as compared to GERG while also being 10 times as fast. The CSP implementation was also more robust, being able to converge successfully more often.


Introduction
Simulations of thermodynamic systems are necessary in the design of new industrial processes. It is therefore essential that the results of such simulations reflect how the system would perform in reality. At the same time, it is an advantage if the thermodynamic simulations do not require much time per simulation, as rapid results enable more options to be investigated in a given amount of time. For optimisation programs it is important that the thermodynamic models are robust and accurate in order for the program to obtain useful information when evaluating the objective function and the constraints. In this study we investigate the trade-offs and compromises between computational time, accuracy and robustness on an organic Rankine cycle (ORC) optimisation problem.
Although thermodynamic models play a major role in simulations of power cycles, little work seems to have been dedicated to investigating how they influence the accuracy, computational time and robustness of the simulation. Because it has been difficult to find such studies for power cycles, we have expanded the focus of our literature survey to include the study of thermodynamic models for other processes as well. However, of all the studies that were found that investigated equations of state, none studied the robustness of the thermodynamic model when implemented in a simulation tool, nor the time it took to obtain results. In fact, they only examined the accuracy of the various thermodynamic models.
Frutiger et al. [1] studied how the simulation result of a power cycle can significantly change as a result of the uncertainty of a few critical parameters such as critical temperature, critical pressure and the acentric factor. Applying a Monte Carlo procedure to 40 different working fluids, they showed how ignoring the uncertainties in the parameters can lead to highly erroneous simulation results, as the 95% confidence interval may be large in cases where important parameters are very uncertain. Furthermore, switching the ranking scheme for the potential fluids between being based on mean net power and according to the lower bound of the 95% confidence interval, changed the relative standing of the fluids. Stijepovic et al. [2] studied how measures of Rankine cycle performance, such as net power, exergy efficiency and cost of power, are influenced by a few parameters of the working fluid, qualitatively showing for example that the molar volume of saturated liquid has a negative impact. Their results highlight the need for accurate calculations on this parameter. Mazzoccoli et al. [3] studied the use of seven different equations of state (EOSs) to predict the temperatures, pressures and densities of binary mixtures with CO2, and compared the results with experimental data. The authors comment that deviations in temperature and pressure from the experimental values are negligible, and so the study focuses on predicting the liquid density. The results show that the cubic EOSs have the highest margin of error, with Lee-Kesler-Plocker (LKP), Perturbed Chain -Statistical Associating Fluid Theory (PC-SAFT) and Groupe Européen de Recherches Gazières (GERG) EoS all offering better predictions. For LKP and PC-SAFT, the authors had to fit the binary interaction parameters, k ij , to obtain the best results. This was unnecessary for GERG, which managed to have among the lowest deviations in this study. The use of different EOSs to predict thermodynamic parameters was also looked into by Li and Yan [4,5], who employed various cubic EOSs to predict the vapour-liquid equilibria (VLE) and molar volumes of liquids and gas for binary mixtures with CO2. Abdollahi-Demneh et al. [6] studied 23 EOSs and evaluated their ability to predict thermodynamic parameters such as vapour pressure and enthalpy of vaporization of 102 pure substances, and found that LKP and the three-parameter cubic EOS from Patel-Teja (PT) were among the best EOSs, having an average absolute percentage deviation of less than 3.41% for the saturated liquid molar volume and less than 1.78% for the saturated vapor molar volume of alkanes and cycloalkanes. Liu et al. [7] studied how well seven EOSs were able to predict the densities of five pure hydrocarbons across a large range of pressures. They also studied implementations of the Soave-Redlich-Kwong (SRK) EOS that used temperature-dependent volume translations, and a volume translation that was both temperature-and density-dependent. This differed from other studies that did not use these dependencies. Their findings show that PC-SAFT consistently has the lowest mean absolute deviation. The study also showed that temperature and density translating the SRK EOS did not much improve the result of the SRK EOS. Orbey and Sandler [8] studied different mixing rules in combination with the Peng-Robinson (PR) EOS for estimating the excess enthalpy and VLE of four binary mixtures. The authors consider this to be a stringent test for the models. They found that fitting the parameters to data for one property, and then predicting the other, yielded poor accuracies, although fitting both simultaneously gave good results for both parameters. Brown [9] studied three different implementations of the PR EOS, and how the accuracy of the EOS changes, using BACKONE and REFPROP 8.0 as references. The three implementations differed in how many parameters were inputs in the calculation, with the remaining parameters having to be estimated. As expected, their results show that the implementation with the fewest estimated parameters gave the most accurate results. While previous studies compared purely simulated results to certain references, Kuboth et al. [10] showed how combining experimental data with the simulations by implementing correlations for heat exchanger performance, pressure loss and refrigerant sub cooling could significantly improve the accuracy of the results gained from the simulations of power cycles.
Other studies have focused primarily on the accuracy with regards to a single or a few parameters of the various EOSs, and while it is very important that the EOSs accurately predict thermodynamic parameters, there are additional characteristics of the EOSs that are important when they are used in simulation or optimization tools. One such aspect is how robust an EOS is in the given simulation tool; if the simulation often crashes or requires frequent restarts, it may be too cumbersome to use in spite of whatever improvements in accuracy it may achieve. Complex models that are highly accurate are often less robust than their simpler alternatives, particularly when calculating phase equilibria. The reason for poor robustness is explained in detail by Wilhelmsen et al. [11]. The crux of the problem is that the functional relationship between two parameters (e.g., pressure and density at a constant temperature) can have several roots in these more advanced models, and therefore there are several infeasible solutions that must be eliminated. Another feature that is important in an EOS implementation is its computational time. Having an EOS that requires less computational time is an advantage, as it allows the user to obtain results rapidly and make adjustments to the model where necessary. These facets have not been discussed and compared in literature to the best knowledge of the authors. Given how little these aspects have been discussed in the literature, we set out to provide quantitative data on all three categories. The EOSs are subsequently ordered based on overall performance. Furthermore, where the published studies have primarily revolved around how accurately EOSs estimate individual parameters, this work focuses on how well the EOSs performed when optimizing an ORC and shows how the use of different models find different sets of optimal operating conditions and different maxima for power production.

Methodology
The EOSs in this work were implemented in an in-house process optimisation framework that maximises the net power output from an ORC. The simple process layout is shown in Figure 1 and consists of a waste heat recovery unit (WHRU), an expander, a desuperheater, a condenser and a working fluid pump. Other details of the process can be found in Table 1. The optimisation method used for solving the problem is a gradient-based, constrained non-linear optimisation routine called NLPQL developed by Schittkowski [12]. More details on how the optimisation framework communicates with NLPQL are described by Skaugen et al. [13]. In configuring the system, we can choose which of the process parameters, hereafter called free variables, that the optimizer may change to maximize the objective function. This study deals with four such variables: The working fluid flow rate, the pump outlet pressure (i.e., evaporation pressure), the expander outlet pressure (i.e., the condensing pressure) and the WHRU outlet temperature. The optimisation routine is allowed to vary these four parameters within defined lower and upper limits. These ranges are shown in Table 2 and are not changed between the individual EOS investigations.  The WHRU was divided into 300 sections of equal enthalpy intervals, and the temperature change across each zone was calculated for the heat source and the working fluid stream in a single pass counter-flow configuration. The local minimum temperature difference between hot and cold streams were returned to the optimisation routine as an inequality constraint, meaning that this must be above a specified minimum value. During the calculation of heat exchanger performance, the minimum temperature difference is allowed to become negative, which will result in the optimisation routine needing to adjust one or more of the process parameters to raise it above the minimum again. This method was also used for the de-superheater and the condenser, but these were divided into 50 sections instead. Our model does not use correlations to estimate detailed heat exchanger performance. Instead, only the thermodynamic state points (specific enthalpy and pressure) at the inlet and outlet of each component were calculated. In the optimisation, the expander operates with a constant isentropic efficiency of 85%, whereas the pump has an isentropic efficiency of 70%. Both of these operate adiabatically between the condenser pressure and WHRU pressure.
There are also several constraints in the model. One of these is that there cannot be two-phase flow through the expander, and this is ensured by requiring single phase vapour at both the inlet and outlet. Furthermore, a minimum temperature difference between the hot and cold stream of 15 K was set in the WHRU, and the heat source was allowed to be cooled to a minimum temperature of 80 • C. This represents a lower limit in a real exhaust gas that may contain moisture and sulphur. The minimum pinch temperature in the condenser is 5 K. Because the model did not constrain the areas of the heat exchangers, the optimal solution always reached the minimum pinch in both heat exchangers.
The EOSs investigated in this work are shown in Table 3. The EOSs were chosen due to their popularity in commercial process simulators and published works. There are a set of two-and three-parameter cubic EOSs using different mixing rules and in some cases using volume shift (Peneloux) for improved density predictions in the liquid/dense phase. Corresponding state principles (CSPs) using different reference equations were also investigated, as were multi-parameter high accuracy EOSs. Details of the different categories of the EOSs, and their implementation, can be reviewed in [14]. Details about the CSP EOSs are explained by Ibrahim et al. [15]. It should be noted that the EOSs prefixed CSP in this work are called SPUNG in [15]. In the notation used here CSP-PR means that the PR cubic EOS is used for scaling temperature and pressure between the actual mixture and the reference fluid. Here propane was used as a reference fluid, because the reference equations-MBWR32 [16] and the newer equation from NIST (C3NIST) [17] -are highly accurate for this fluid. Table 3. Equations of state investigated in this study.

Cubic Equations of State Multiparameter Model
Peng-Robinson (PR) [18] GERG [19] Soave-Redlich-Kwong (SRK) [20] Corresponding state principles (CSP) PR-Peneloux [21] Lee-Kesler (LK) [22] PR-UNIFAC [23] CSP-PR [15] with MBWR32 [16] as reference EOS PR-UMR (Universal Mixing Rule) [24] CSP-PR [15] with C3NIST [17] as reference EOS PR-HV1 [25] Statistical Associating Fluid Theory SRK-Peneloux [26] PC-SAFT [27] Patel-Teja (PT) [28] The objective function for the optimisation is the net power from the cycle, where the subscripts refer to the various state points in Figure 1, with h being the specific enthalpy and . m w f being the working fluid mass flow rate. The optimisation model uses numerical two-sided differentials when calculating gradient information for NLPQL as it seeks the optimum values for the free process parameters. This method is highly dependent on the accuracy and robustness of the underlying thermodynamic model chosen, as models that are not robust may lead to the optimiser giving very different results between the different optimisations.
In evaluating the accuracy and time of the optimisation, all the equations of state are set to optimise an ORC under the same initial conditions. GERG is used as a reference for the accuracy of the EOSs, as this equation predicts the properties within the experimental accuracy for mixtures of propane and butane [29]. The computational time is the measured time taken until the optimisation is finished, as reported by the program itself. It is important to note that all the EOSs are implemented identically in this model, ensuring that any differences in time are a result of the EOSs themselves, and are not due to differences in implementation. The model is run on a single processor core.
For the analysis of the robustness of each EOS, the optimisation problem was repeated 50 times, starting from random initial conditions (within the specified permissible range) for the four free process variables. The random initial conditions were generated using the same seed, ensuring that all the EOSs had the same initial data set. Performing multiple optimisation runs also increases the probability of finding the global maximum for the problem. The maximum recorded net power obtained from all of the 50 optimisations was marked as best, and to evaluate the consistency of the other 49 results, the results were given a score depending on how close they were to the best result. This score system is detailed in Table 4. Each successful optimisation gave 10 points, rewarding EOSs that consistently gave a result and penalising EOSs that tended to crash. This means that the highest score for an EOS is 843; 500 for 50 successful optimizations, and 7 points to each optimization that is very close to the optimal result, totalling 343.

Results
The results of the analysis are given in terms of their accuracy, computational time and robustness.

Accuracy
Each EOS was used in the optimisation program to find the set of variable parameters that maximized the net power from the ORC. This result was compared with the optimised result found with the GERG EOS. The percentage deviation of optimized performance parameter X was found using where X is a generic parameter. The deviation in maximum net power for each thermodynamic model used in this work is summarized in Figure 2. The differences between the calculated maximum net power obtained may have more than one explanation. One of these is that the optimizer finds different optimal state points when different EOSs are used. This is evident when investigating the temperature at the inlet of the expander for each EOS, as is shown in Figure 3, and the working fluid flowrate, shown in Figure 4. Here, the cubic EOS consistently report temperatures that are around 8 K higher than the reference while the flow rate is about 5-6% lower. There is no improvement including the volume shift or using a 3-parameter cubic equation like the PT EOS. Switching to corresponding state methods or PC-SAFT improves the accuracy significantly.  Given that the optimized process parameters play a major role in the optimized solution, it is important also to include the deviations from the optimal process parameters reported by the reference, rather than just the net power. To do this, we define the total deviation, which includes all the free process parameters, where P 3 is the expander inlet pressure, T 3 is the expander inlet temperature, T 5 is the condensation temperature, and . m w f is the working fluid mass flow rate. Figure 5 shows the calculated total deviation of each EOS.  Figure 5 shows that combined for all the free process variables, all the two-parameter cubic EOSs show a total deviation from the reference of between 8% and 11%, with no improvement when volume shift for density calculations is included. The three-parameter cubic EOS from PT does not offer a significant improvement in the total deviation either. With the CSP methods the result improves significantly, reducing the total deviation to less than 1%.

Computational Time
As mentioned above, how quickly the different EOSs reach a solution is an important consideration when choosing an EOS. Figure 6 shows the relative computational times of the individual EOSs, with GERG as the reference, i.e., GERG has a relative computational time of 100. These were calculated using GERG is also chosen as a reference in this analysis for time because this makes it easier to investigate the trade-off between computational time and accuracy or robustness.

Robustness
While accuracy and computational time are both important, it is also highly beneficial if the EOS is robust and it is essential that it is consistent. What is meant by "robust" is that the EOS does not crash or otherwise require restarts when searching for the optimized solution, and having a consistent EOS means that it finds the same solution regardless of initial conditions. This analysis uses the point scheme detailed in Table 4 for 50 optimisation runs using random initial conditions. Figure 7 shows how well each of the EOSs performed. Surprisingly, GERG scored relatively poorly on both robustness and consistency in this investigation, being bested by several of the simpler EOSs. Not all of the optimisation runs completed successfully, and there was a significantly larger spread in the results among those that did. The latter may be a result of more relaxed tolerances in the implementation of the GERG EOS in REFPROP 9 than in our in-house implementation of the other models. The LK EOS is a two-fluid corresponding state method and turned out not to be robust for this system with evaporation relatively close to the critical pressure for the propane-butane mixture.

Overall
With the results for each of the three categories in hand, we can order the EOSs to find which seem to have a good overall performance with regards to accuracy, computational time and robustness. To order the EOSs, first the performances of each EOSs in every category were normalised to give the best performer a score of 1, and the worst 0. The other EOSs are distributed according to how they performed compared to the best EOS in that category. An EOS that scores best in each category would thus get the top score of three points. Figure 8 shows the overall performance of each EOS. It is evident that while GERG is accurate, it is heavily penalized for its relatively high computational time. On the other hand, many common cubic EOSs such as PR and SRK display poor accuracy, while being consistent and rapid. PC-SAFT, an EOS that has been identified as accurate, does not seem particularly accurate in this study of a system with light hydrocarbons. This may be a result of not having tuned the adjustable parameters of the EOS. It is ranked the poorest of the EOSs in this work because of its lack of accuracy, in conjunction with being the slowest of all the EOSs. The EOSs that we ranked highest are the two CSP-cubic EOSs. They offer close to the accuracy of GERG but are much faster while also being robust.

Conclusions
We studied the use of several different EOSs in an optimization problem for organic Rankine cycles, and compared the EOSs on accuracy, time and robustness. It was found that for the purpose of simple cycle modelling and optimisation of the scenario in this work, thermodynamic models based on the corresponding state principle (CSP) using a highly accurate reference equation and a simple 2-parameter cubic EOS for the scaling of mixture properties to the reference property offered a good trade-off: Retaining high accuracy while being both fast and robust. The CSP models with propane as the reference fluid were able to reach an optimal solution roughly 10 times faster compared to the GERG EOS, while only deviating by 0.7%. Moreover, the CSP models did not fail any of their optimisations, whereas GERG failed a few. The performance of the CSP EOSs with regards to computational time, robustness and accuracy culminate in them seeming to be well-suited to the simulation and optimisation of simply modelled power cycles using hydrocarbon mixtures as the working fluid and operating in the vicinity of the critical pressure, but a similar study using several working fluids must be performed in order for this to be confirmed. Additionally, these results were found using only one set of parameters for the heat source and heat sink. Before they can be accepted as general to all organic Rankine cycles, it must be investigated if other sets of such parameters attain the same results.
Author Contributions: G.S. supervised the work and developed the model used; G.S. also acquired the data used for the analyses in the work; G.D. analysed the data, created the figures and wrote the original draft; and both G.S. and G.D. reviewed the manuscript internally, whereupon G.D. modified the manuscript.

Funding:
The authors gratefully acknowledge the financial support from the Research Council of Norway (EnergiX program, grant no. 255016/E20) for the COPRO project, and the user partners Equinor, Hydro, Alcoa, GE Power Norway and FrioNordica.

Acknowledgments:
In this section you can acknowledge any support given which is not covered by the author contribution or funding sections. This may include administrative and technical support, or donations in kind (e.g., materials used for experiments).