This section describes in detail the overall architecture of the simulation model as well as each independent component model of which the cycle model is composed. The most general case of an ORC with internal heat exchanger is considered in this section.
2.1. Thermo-Physical Property Wrappers
Most of the computational effort of a detailed cycle model takes place while retrieving and evaluating the properties of the working fluids. Therefore, the property wrapper must combine computational efficiency with accuracy and a sufficiently broad database of fluids.
Both CoolProp 5.1.2 [
26] and REFPROP 9.1 [
27] have been included within the detailed ORC model to cover the most common ORC applications in terms of refrigerant working fluids and their mixtures, as well as hot and cold source fluids, compressible and incompressible. In particular, incompressible fluids such as thermal oils or coolants are only available in CoolProp. A modified version of the
Refprop2.py [
28] has been implemented which allows us to call the properties from REFPROP with the same syntax as in CoolProp. For consistency with CoolProp, the
PropsSI function has been used, with only SI units. The possibility to call refrigerant mixture properties based on mass fraction has also been added to the original version of
Refprop2.py.
Positive displacement expanders require lubrication to ensure the correct running of the rotating parts and to reduce the internal leakage paths. To the authors’ knowledge, lubricant oil property libraries are not readily available open-source. An additional property wrapper has been added to calculate both the properties of the lubricant oil and the mixture of oil and refrigerant. Many of the oil properties were obtained from [
29]. To simplify the analysis, the mixture has been treated as a pseudo-pure fluid,
i.e., an ideal mixture model has been used. For a homogeneous mixture of two components with equal velocities, mixture properties such as specific enthalpy, specific entropy, specific internal energy and exergy, the specific volume, and the mixture constant pressure specific heat can be calculated by introducing an oil mass fraction,
. For example the mixture specific enthalpy is given by:
The other aforementioned properties can be calculated in a similar fashion. The mixture density and mixture thermal conductivity are calculated based on a void fraction weighted average correlation:
where the homogeneous void fraction,
α is given by:
with the slip ratio
S equal to unity. Many correlations exist to calculate the mixture viscosity. Two models have been implemented. The first one is based on the model proposed by McAdams [
30]:
The second one is the correlation developed by Dukler
et al. [
31] based on the mass averaged value of the kinematic viscosity:
The solubility of refrigerants in lubricant oils is also included and is discussed further in
Section 2.7.
2.2. Plate Heat Exchanger Model
The heat exchangers of the overall ORC model (typically evaporator, condenser and regenerator) are modeled with a general steady-state moving boundary approach. In this work only plate heat exchangers (PHEX) are considered. However, the model can be readily adapted to other types of heat exchangers by modifying the geometric parameters and heat transfer and pressure drop correlations. The model is nearly identical to the one proposed by Bell
et al. [
32] and available open-source within the ACHP project, version 1.5.0 [
33]. In the following the general algorithm of the model is briefly described and the modifications highlighted.
The modeling scheme considers a steady-state, counter-flow, moving boundary method which has been generalized to account for any phase condition for either side, namely hot side and cold side. The main inputs of the model are the PHEX geometry (number of plates, dimensions of the plate, wave length, amplitude of corrugation and Chevron angle, thermal conductivity of the plate) and the mass flow rate and inlet conditions of both streams in terms of pressure and enthalpy. A degree of freedom consists of the possibility to assign an extra channel to the hot or cold stream in the case of an even number of plates,
i.e., odd number of channels. In fact, in a PHEX with
N plates, there are
total channels, which results in
active counterflow sections, as shown in
Figure 2a. The details of the calculation of the PHEX surface area can be found in [
33]. The moving boundary heat exchanger model is characterized by its use of separate zones, the boundaries of which are defined by the thermodynamic phase change points of each fluid. The most general case where both fluids enter with single-phase states, undergo a complete phase change, and exit with single-phase states is shown schematically in
Figure 2. The heat exchange rate with the surroundings is neglected in the current version of the model.
The physical bounds of the total heat transfer rate are obtained based on a pinch point analysis. Only phase-change pinch points are considered in the moving boundary model and only sub-critical conditions are allowed. The maximum physical heat transfer rate,
, is first assumed to be given by:
Once
is evaluated, a list of enthalpies at the inlet, exit, and zone boundaries (if any) of each working fluid is constructed. Where a phase boundary is found for one fluid, the corresponding state in the other fluid is determined using an energy balance. Pinch point violations are checked and handled, as described in [
32].
For a given heat transfer rate and inputs, the number of zones and phase boundaries can be directly determined. The ability of the code to correctly handle any phase configuration is useful for ORC systems. An extension to the code proposed in [
32] has been added to handle incompressible fluids, such as thermal oils and lubricant oils where only the liquid phase is allowed.
For each zone, the surface area required to achieve the heat transfer rate defined by the zone boundaries is determined. The required heat conductance for the zone, , is calculated by applying the log mean temperature difference (LMTD) method. The use of the LMTD method, rather than the effectiveness-NTU method, made the code more concise and still able to handle two-phase conditions for pure fluids or mixtures.
The required heat conductance for a zone is defined by:
Zero division errors that may occur while evaluating the LMTD function in Equation (
8) have been avoided by introducing numerical tolerances. The heat transfer area fraction of the heat exchanger associated with the required heat conductance can be calculated. Noting that for heat exchangers with constant cross-section along their length, such as PHEX, the heat transfer area fraction reduces to a length fraction defined by:
where
is the center-to-center distance between the ports of a plate and
P is a constant dimension of the surface perpendicular to the direction of the heat transfer. A bounded solver is used to iterate on the heat transfer rate until the sum of length fractions in each zone is equal to unity, indicating that the entire surface area of the heat exchanger was utilized. In particular, the convergence is attained by driving to zero the following residual function:
The bounded solver used in the current version of the model is based on Brent’s method. The brentq function available in the scipy.optimize module for the Python programming language has been employed in the code.
Due to the fact that the moving boundary model presented is general, the same code can be used to simulate the evaporator, condenser and regenerator. Proper heat transfer correlations for evaporating and condensing streams have been implemented. For instance, the Cooper pool boiling correlation is used to find the heat transfer coefficient for evaporating flow [
35]. To allow for the possibility of evaporating flow at high quality, the Shah chart correlation for boiling heat transfer has also been included and it can be chosen optionally by the user [
36]. The heat transfer coefficient correlation for condensing flow used in the model refers to the study conducted by Longo on condensing refrigerant flow in PHEX [
37]. Note that no general two-phase heat transfer correlations have been implemented for two-component mixtures. The pressure drops of each fluid are calculated by assuming an appropriate pressure drop correlation for each zone. In the solution scheme, the pressure drops are not included in the heat transfer calculations but they are calculated independently. This makes the model solution more straightforward by eliminating the need to iterate between pressure drop and heat transfer in each zone. Additional documentation can be found in [
33].
The PHEX moving boundary model allows the estimation of the refrigerant charge in the heat exchanger. Due to the fact that during the operation of an ORC, charge fluctuations occur in the PHEX, a detailed model to estimate the charge in each zone has been developed. In particular, the total working fluid charge in a PHEX is given by:
where
is the combined volume of all channels for either the primary or secondary fluid side,
represents the volume of a single zone and
is the average fluid density of the zone. In the case of single-phase zone,
is calculated assuming an average temperature between inlet and outlet of the zone. For two-phase zones,
is obtained by weighting the density of the liquid phase and vapor phase with an average void fraction for the considered zone. Mathematically, this can be expressed as:
where the average void fraction over the length of the zone is given by:
where
is the slip model defined by Zivi [
38]. The main inaccuracies of this charge model may lie on the void fraction model used to estimate the density in the two-phase zone, where inaccuracies may occur. An improvement to this model proposed in this paper would be to implement different void fraction models to account for different flow patterns [
38]. Another source of error is that a quality-averaged void fraction is used rather than a length-averaged void fraction. This is done for convenience to avoid estimating the quality variation along the length of the heat exchanger, but it does not necessary represent the distribution of voids within the zone.
Examples of temperature profiles across the PHEX (both as condenser and as evaporator) are shown in
Figure 3. The fluid combinations are representative of the two ORC installations considered in
Section 4. The PHEX model allows us to make use of the extended incompressible fluid library available in CoolProp [
26]. Three examples are proposed: a condenser with a mixture of ethylene glycol and water as cooling medium,
Figure 3b, an evaporator using a thermal oil (Therminol 66) as hot source,
Figure 3c, a lubricant oil heater,
Figure 3e. An example of heat transfer between a zeotropic mixture of R134a and R245fa and hot water is presented in
Figure 4a. Phase-change conditions are also allowed in both streams. In ORC systems, this situation may occur at the regenerator, as shown in
Figure 4b.
2.3. Pump Models
The pump models implemented are based on empirical maps. A distinction is made between volumetric and centrifugal pumps. Due to the fact that volumetric and centrifugal pumps have different characteristic curves, appropriate empirical correlations are selected to better capture their behaviors. In both cases, two empirical correlations have been used. One models either the volumetric or mass flow rate, and another correlation gives the performance of the pump, either the input power or the pump efficiency.
For the volumetric pump two sets of correlations are considered:
- (i)
mass flow rate and input electric power;
- (ii)
normalized volumetric flow rate and pump efficiency.
In the case of the centrifugal pump, only one model has been implemented:
- (i)
mass flow rate and pump efficiency with normalized input values.
The different models for each pump type are listed for compactness and clearness in
Table 1. The typical volumetric pump curve is characterized by a linear relationship between the volumetric flow rate,
and the rotational speed,
for different discharge pressures,
. A linear correlation can be obtained to express the volumetric flow rate of the pump as a function of the rotational speed and discharge pressure:
where
,
,
and
are known from the pump manufacturer curves. From Equation (
14), simplified empirical correlation can be derived and coefficients can be fit to available experimental data. By considering the mass flow rate rather than the volumetric flow rate, Equation (
14) can be rewritten as:
where the parameters
a and
b are obtained by minimizing an objective function. By introducing a set of fitting parameters,
,
,
and
, and by normalizing the input and output values with respect to a set of reference values to improve the robustness of the regression, an equivalent form of Equation (
14) can be obtained:
where the discharge pressure
in Equation (
14) has been substituted by the pressure difference across the pump,
, due to a better agreement with the experimental data. The centrifugal pump curves are well represented by second-order polynomial functions. By choosing the pressure increase provided by the pump,
and its frequency,
, as input variables, the volumetric flow rate can be expressed as,
where the coefficients
are obtained by regression.
The input power of the pump,
, is calculated as the ratio of the reversible work rate,
, to the pump efficiency,
. A choice can be made whether using an empirical correlation for the input power or the pump efficiency and to calculate the other term. In the case of the volumetric pump both correlations have been developed. By defining the reference power of the pump,
, as a linear function of the pump pressure difference at the chosen reference rotational speed,
, the pump power at any rotational speed can be calculated as:
A correlation for the pump isentropic efficiency can be used instead of the pump power because it provides better physical constraints on the model and it guarantees that the Second Law of Thermodynamics is not violated. Analogously to the pump power equation, the pump efficiency can be related to a reference pump efficiency,
at a chosen rotational speed. From experimental data, an exponential decay of the pump efficiency with the pump speed was observed. The functional form is given by
where the reference pump efficiency is a function of the mass-specific reversible work of the pump,
at the reference pump speed, and the coefficients
are determined empirically. Mathematically:
By introducing Equation (20) in Equation (19) and by normalizing the input values with respect to the reference conditions, the final form of the correlation for the pump efficiency can be obtained, as listed in
Table 1. A similar approach has been considered for the centrifugal pump. The pump efficiency has been correlated to different non-dimensional terms normalized with respect reference values. Due to the fact that experimental data of the pump isentropic efficiency did not show a clear trend with respect to the mass flow rate, rotational speed and pump pressure difference, a third-order polynomial function with several cross-terms has been selected. The pump isentropic efficiency has been expressed in terms of the pressure difference, pressure ratio, density of the working fluid and frequency. The functional form can be written as:
where
and
. The complete formulation can be found in
Table 1.
2.4. Expander Models
In this study, positive displacement expanders are considered as suitable devices for low-grade heat recovery [
6]. The expander has been modeled using an empirical map and also by means of a semi-empirical model.
Empirical model. The performance of a volumetric expander is well characterized by the filling factor and the isentropic efficiency. Two empirical correlations can be developed to account for different operating conditions of the expander. The filling factor along with the displacement rate of the machine is used to calculate the mass flow rate for a given inlet state conditions and working fluid. It is defined as:
where
is the displacement volume of the expander. An empirical correlation similar to the one proposed by Declaye
et al. [
39] is employed. Its expression is given as:
where
are the fitting parameters to be obtained by regression. Note that the filling factor is a function of the inlet pressure, the pressure ratio and the rotational speed which have been normalized by a set of reference values as follows:
A second empirical relationship is used to model the expander isentropic efficiency. By doing so, the expander power output can be obtained by calculating the mass flow rate from Equation (
22) and the mass-specific isentropic work for a given expender inlet temperature and pressure ratio. As first proposed by Declaye
et al. [
39], Pacejka’s equation is suitable to represent the behavior of a volumetric expander isentropic efficiency with respect to either the volume ratio or pressure ratio at different rotational speeds. Pacejka’s equation has been applied to scroll expanders [
39,
40] as well as to a single-screw expander [
22] with high accuracy in the prediction of the isentropic efficiency. This shows the versatility of the equation as a modeling tool for different types of volumetric expanders. The original equation described by Pacejka and Bakker [
41] is expressed by:
where
y and
x indicate a position on the ordinate and abscissa, respectively. Physical values related to the expander operating conditions are selected for the abscissa and ordinate axes to represent the expander isentropic efficiency. Furthermore, functional forms are required for the coefficients appearing in Equation (
25). A more in-depth description of the equation and its parameters can be found in [
34,
39].
Semi-empirical model. A semi-empirical model has a computational time comparable to an empirical model while retaining physically meaningful parameters that can be identified from experimental data. The typical increase in robustness due to less modeling detail involved and the absence of differential equations makes the semi-empirical model suitable to be included in an overall cycle model. The semi-empirical model described here is an extension to the well-known semi-empirical model initially developed by Winandy
et al. [
42] and then successfully applied to scroll expanders by Lemort
et al. [
43] and Quoilin
et al. [
1]. Such a model has been demonstrated to be applicable to different types of volumetric expanders, e.g., scroll, screw and piston types, by Guillaume
et al. [
9]. The expander model is modified to take into account the presence of oil in the machine throughout the expansion process. In the following, the major changes related to two-phase pressure drops, leakage and the actual expansion compared to the original model are only briefly described. A more comprehensive description can be found in [
40,
44].
The evolution of the working fluid and lubricant oil through the expander is decomposed in several steps, as shown in
Figure 5:
adiabatic supply pressure drop;
isobaric supply cooling;
fictitious internal leakage;
adiabatic reversible expansion to the internal pressure ratio imposed by the built-in volume ratio of the expander;
adiabatic expansion at a constant volume;
adiabatic mixing between the flow exiting the expansion chamber and the internal leakage flow;
isobaric exhaust heat transfer;
adiabatic discharge pressure drop.
The working fluid and lubricant oil are considered to be in thermal and mechanical equilibrium so that a homogeneous mixture model can be applied. The total mass flow rate entering the expander is given by
where
represents the ratio between the oil mass flow rate to the working fluid mass flow rate. The pressure drop through the suction port of the expander is modeled by calculating the mass flow rate of a compressible gas-liquid mixture through a nozzle. Based on the Chisholm [
38] model, such a mass flow rate can be calculated by imposing the conservation of momentum between upstream and downstream conditions of the orifice as:
where
represents the two-phase flow coefficient calculated according to Morris [
45], the effective specific volume or momentum specific volume
is calculated by assuming the case of a separated-phase flow with some liquid entrainment in the gas phase and flowing at the same velocity. The correlation of
and the associated effective slip ratio between the phases can be found in [
38,
45]. The coefficient
is the throat to upstream area ratio of a flow path, which can be approximated to be zero if the throat area is treated as being infinitely smaller than of the upstream area. Under this assumption, the lowest flow rate possible with this model is computed [
29].
Dealing with two-phase leakage flows leads to complex fluid dynamic effects. In order to take into account the presence of oil in the leakage paths without increasing the computational effort of the model, the total leakage area has been corrected by a factor proportional to the oil ratio that enters the machine. In other words, the more the oil fraction is increased, the more the working fluid leakage is reduced, which is indeed one of the purposes of the lubricant oil. Mathematically, this statement is expressed by still considering a compressible flow of working fluid through a nozzle where the throat area,
, is a parameter to be determined empirically:
The working fluid leakage flow rate is given by:
Therefore, the internal refrigerant mass flow rate that contributes to generate positive work is given by:
and the new oil fraction entering the expansion chamber is defined as:
The expansion process is divided into two steps as shown in the original model [
1,
43]. At first, the oil-refrigerant mixture undergoes an isentropic expansion where the specific entropy of the mixture defined as
is kept constant. The internal isentropic specific-work is given by:
Due to the presence of oil, the built-in volume ratio becomes an effective built-in volume ratio:
The second step is an expansion at constant machine volume. The work associated with this second step is expressed by:
where implicitly it is assumed that
. Finally, the overall isentropic efficiency of the expander is given by:
where
accounts for the visco-mechanical losses in the expander. The visco-mechanical losses are assumed to have a linear dependency with the rotational speed and the oil fraction. A constant frictional torque,
, is introduced in the model along with a viscous torque due to the presence of lubricant oil. The viscous torque is given by:
where
represents the viscous coefficient. Hence, the total losses are expressed as
The semi-empirical model is closed by a shell energy balance over a fictitious wall envelope:
The semi-empirical model includes a set of fourteen parameters that have to be identified. A multi-parameter optimization is carried out by minimizing an objective function based on the errors on working fluid mass flow rate, the discharge temperature and the expander power output. These represents also the three main outputs of the model along with the expander isentropic efficiency. The objective function is defined as a weighted sum of the relative error between predicted and measured values for each of the variables considered:
The weighting factors are chosen to give each error term approximately the same weight. For the expander discharge temperature term, the weight factor can be expressed as the difference between the maximum and minimum discharge temperature values. However, in this work
has been selected for both expanders. A similar objective function was also used by Quoilin
et al. [
1].
2.6. Line Set Model
The line sets of the ORC system carrying the working fluid are considered in the overall cycle model. The main reasons to consider the line sets are to determine the pressure drop between the main cycle components and to give a more accurate estimation of the charge of the working fluid. The following line sets are included in cycle model:
liquid receiver exit to pump inlet;
pump exit to liquid side of the regenerator;
regenerator liquid side exit to evaporator inlet;
evaporator exit to expander inlet;
expander exit to vapor side of the regenerator;
regenerator vapor side exit to condenser inlet;
condenser exit to liquid receiver inlet.
Each line set is modeled as an equivalent tube with inner and outer diameters (ID and OD). The equivalent length of each section (
) of the line set includes straight sections, tees, elbows and other fittings. It is calculated using the method and loss coefficients of Munson
et al. [
46]. In the calculation of the pressure drops, a distinction is made to account for single and two-phase flow conditions, as proposed by Woodland [
34]. In the case of single-phase flow, the Darcy friction factor is used and the pressure drop for a certain ID is given by:
If two-phase flow conditions occur, the original form of the Lockhart-Martinelli method for two-phase frictional pressure drop in tubes [
47] is used. The equivalent lengths of the fittings are computed using Churchill’s correlation for the friction factor under the hypothesis of liquid-only. The pressure drop of each section of the lineset is then given by the product of the Lockhart-Martinelli pressure gradient and the total equivalent length of the lineset:
The total pressure drop of a certain lineset is given by the sum of each section’s pressure drop.
The lineset model includes also a simplified model for the heat losses. The overall heat transfer coefficient of the lineset,
, is the sum of the thermal resistance of the tube, the thermal resistance associated with the insulation and the convective
values associated with the inside and outside of the tube:
If heat losses exist in a lineset with two-phase flow, then the pressure gradient in Equation (
42) represents an average pressure gradient along the lineset flow direction.
The working fluid charge of the lineset is calculated analogously to the method applied to each zone of the heat exchanger. If a lineset is assumed to be adiabatic, then the average void fraction (Equation (
13)) under two-phase flow conditions is equal to the void fraction at the inlet quality of the lineset.