Validation and analysis of organic Rankine cycle dynamic model using zeotropic mixture

Energy 197 (2020) 117003

Contents lists available at ScienceDirect

Energy journal homepage: www.elsevier.com/locate/energy

Validation and analysis of organic Rankine cycle dynamic model using zeotropic mixture Jinwen Cai, Gequn Shu, Hua Tian*, Xuan Wang**, Rui Wang, Xiaolei Shi State Key Laboratory of Enigne, Tianjin University, No.92 Weijin Road, Nankai District, Tianjin, 300072, China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 19 July 2019 Received in revised form 22 November 2019 Accepted 18 January 2020 Available online 5 February 2020

Organic Rankine cycle (ORC) is being widely researched in the application of waste heat recovery and renewable energy utilization. The ORC using zeotropic mixtures (MORC) attracts much attention owing to the ability to match better in heat transfer. The system with fluctuant heat source often encounters challenges of safety and efficiency and this issue drives the development of dynamic research. Modeling methods including finite volume (FV) and moving boundary (MB) can well predict the dynamic behavior of system with pure fluids, but the situation with zeotropic mixtures are unclear because of composition shift effect. In this case, this paper establishes the dynamic model of MORC by these two methods separately and compares the simulation with experiment. Results indicate that both methods have similar high accuracy and the composition shift effect in MORC dynamic simulation could be acceptably ignored. The following study on dynamic characteristics shows the response speed of ORC system with R134a/R245fa tends to increase with mass fraction of R134a. The system response speed with the disturbance of fluid mass flowrate is the fastest comparing to exhaust temperature and flowrate. This work contributes to provide abundant reference for dynamic modeling of MORC and lay a foundation for the subsequent control design. © 2020 Elsevier Ltd. All rights reserved.

Keywords: Organic Rankine cycle Finite volume method Moving boundary method Zeotropic mixture Experimental validation Dynamic response

1. Introduction The energy crisis and environmental pollution has brought challenges to the sustainable development of society. Driven by the polices of energy conservation and emission reduction, many technologies have been proposed to meet increasingly stringent regulations, such as waste heat recovery and utilizing renewable sources. In literature research, organic Rankine cycle (ORC) is suitably and widely applied to recover waste heat from industry [1] and internal combustion engine [2,3]and utilize solar [4] and geothermal [5] energy. The ORC technology has attracted so much attention due to its mature theory, high reliability and strong adaptability. According to Montreal protocol [6], traditional pure refrigerants originally adopted in ORC are gradually banned or restricted for the effect of ozone depletion and global warming. The remaining options are hydrofluoroolefins (HFOs) and hydrocarbons (HCs), while

* Corresponding author. ** Corresponding author. E-mail addresses: [email protected] (X. Wang).

(H.

Tian),

https://doi.org/10.1016/j.energy.2020.117003 0360-5442/© 2020 Elsevier Ltd. All rights reserved.

[email protected]

they are flammable. In this case, zeotropic mixtures [7e10], formed by different pure fluids, are introduced to achieve the improvement of system safety [11] and environment performance [12]. Moreover, mixtures are promising substitutes comparing with present pure refrigerants due to non-isothermal phase transition and flexible selectivity. The non-isothermal properties during evaporation and condensation are helpful for better heat transfer match with heat source and hence better system performance. The flexibility refers to that mixtures could be formed from different components with different fraction. Abadi et al. [13] tested the system performance of ORC with R245fa/R134a (0.6/0.4, mole concentration) and pure R245fa, respectively. The experimental results showed that the output power of mixture is 0.3 kW higher than that of pure R245fa with 90
C heat source. Shu et al. [14] investigated the performance of trans-critical Rankine cycle system for waste heat recovery of diesel engine. They selected CO2 mixtures and pure CO2 as working fluids and concluded that net output power of system with CO2/R32 (0.7/ 0.3) increases by 8.8% comparing to that with pure CO2. Wang et al. [15] carried out an experimental research on low-temperature solar Rankine cycle system and selected three fluids: R245fa, R245fa/ R152a (0.9/0.1) and R245fa/R152a (0.7/0.3). Their work resulted

2

J. Cai et al. / Energy 197 (2020) 117003

that average overall efficiency of M1, M2 and M3 is 0.88%, 0.92% and 1.28%, respectively. Fluctuations can usually be observed when ORC is applied to absorb heat from exhaust and solar energy. According to Shu’s research [16], the parameters of exhaust fluctuate over a larger range when the gasoline engine is operated in Highway Fuel Economy Test Cycle. The temperature varies from 400 to 700 K and flowrate changes from 0.01 to 0.04 kg/s. Wang et al. [15] recorded the global solar radiation from 8:00 to 16:00 in spring and found it varies from 500 to 1000 W/m2. The fluctuation of heat source can strongly influence the system efficiency and safety. Liu [17] adopted off-design model to research the engine load effect on ORC performance and found the maximum net power output decreases from 87.8 kW to 22.1 kW when the engine load varies from 100% to 40%. Benato et al. [18] investigated the ORC combined with gas turbine and found severe load changes (0.4e1.0 MW s1) can cause fluid temperature exceed its decomposition limit. These challenges drive the demand of reasonable and effective controller [19,20] for maintaining efficient and stable operation and the control design are generally developed based on the system dynamic behavior. The dynamic model is preferred to simulate the system transient behavior. Current system dynamic model mainly focuses on heat exchanger while considering pump and turbine as static for their small capacity and rapid revolution speed. Hence the dynamic ORC model could be classified into two kinds: model based on finite volume method (FV-M) [21e23] and model based on moving boundary method (MB-M) [20,24e26]. The former divides the heat exchanger into multiple cells along the flow direction and calculates the fluid state by solving the mass and energy conservation equation, while the latter tracks the length of the phase regions such as subcooled, two-phase and superheated. Dynamic ORC models using pure fluids have been well researched in literature and can be useful reference to the modeling of MORC. Wei et al. [27] verified the system dynamic model based on FV-M and MB-M by experimental data and the selected working fluid is R245fa. They concluded that both models can accurately capture the leading dynamic behavior of system, FV-M is more accurate and MB-M is faster for computation and more suitable for control design. Vaupel et al. [28] assessed the two approaches in modeling heat exchangers and selected the measurement data of ORC with ethanol for validation. They found that the FV-M with sufficient cells results in a similar error as the MB model. Shu et al. [22] developed an MB-M of ORC recovering exhaust heat from natural gas engines. They compared the response speed of 14 kinds of pure fluids under the disturbance of fluid mass flowrate and exhaust temperature. Their work indicated that system generally responds slower with the higher critical temperature. Quoilin et al. [16] established a FV-M for small-scale ORC using R245fa. They proposed three different control strategies when the temperature and flowrate varies and concluded that the model predictive control strategy based on the steady state optimization shows the highest overall efficiency (6.6%). Actually, there are some differences between zeotropic mixtures and pure fluids when adopted in ORC system. The phenomenon of composition shift can be observed during the evaporation and condensation process. This effect could cause actual circulation concentration deviation from the charge concentration and low output power [29]. Zhou’s research [30] on R227ea/R245fa shows that component proportion varies from initial 50.60%/49.40% to operational 52.00%/48.00%. Moreover, the temperature of mixtures during phase transition is non-isothermal. The temperature glide leads to greater mass transfer resistance and loss of effective overheating, and therefore weakens evaporation [31]. These complex physical mechanisms increase prediction difficulty of heat transfer correlations. Zhang’s [31] study suggested the improved

correlation still shows a mean deviation of 24.6% comparing to 2091 experimental data points. Considering these differences, it is necessary to build dynamic model for MORC. Current models dealing with composition shift require multiple iterative calculation at steady state and this means coupling composition shift effect into the system dynamic model is not practical. Moreover, study about composition shift occurring during the heat transfer process of MORC reveals that the concentration variation is relatively small (<4%) [29e31]. Therefore, this paper attempts to firstly establish MORC dynamic model with ignoring composition shift, then evaluate the reasonability of this simplification based on experiment results. Since zeotropic mixtures are different from pure fluids and the dynamic process of ORC using mixtures is rarely researched in the literature, this paper concentrates on the dynamic modeling of MORC. Considering the finite volume method and the moving boundary method are often used in pure fluid dynamic modeling, this paper establishes two system dynamic models using those methods separately in order to figure out whether they can also be applicable to mixtures. Subsequently, a kW-scale ORC test bench is introduced to carry out the dynamic experiment. Experimental verification shows that two methods are suitable for MORC dynamic modeling. The high smulation accuracy indicates the acceptance of ignoring composition shift in MORC dynamic modeling. The comparison between FV-M and MB-M indicates that the former is more applicable due to the capacity of dealing with phase vanishing and forming. The remain of the paper is organized as follows: Firstly, both FVM and MB-M for MORC are presented in section 2. Moreover, the problems encountered in modeling are discussed. In section 3, a kW-scale ORC test bench is introduced and the uncertainty has been calculated for ensuring the reliability of experiment. Section 4 provides the validation of these two system models based on experiment. After that, the FV-M is used to analyze the response characteristics of system with R134a/R245fa in section 5. The conclusion is attached in section 6. This work could contribute to provide reference for the dynamic modeling of MORC and lay a foundation for the subsequent system control design. 2. Model description Fig. 1 illustrates the conceptual scheme of the ORC system. The idea of establishing system-level model is deriving sub-component models firstly and then combining them based on the sequential relationship. The heat exchanger model, as a key component in system [32], is introduced firstly. Remaining models of pump, turbine and liquid receiver are then presented. The model construction is implemented in MATLAB [33] and

Fig. 1. Scheme diagram of the ORC system.

J. Cai et al. / Energy 197 (2020) 117003

SIMULINK. Adopting the built-in S-function to generate subcomponent models and then, using the blocks to form system model. Meanwhile, the required thermodynamic properties of mixture during the simulation process can be calculated by calling REFPROP v9.0 [34]. Assumptions are indispensable before deriving model equations and they can be found in literatures [20,25,35]. 1) The heat release to the surrounding is ignored. 2) The momentum conservation is omitted in models and therefore, the pressure is spatially uniform over the heat exchanger. 3) The flow in exchanger is one-dimensional and horizontal. 4) All the heat transfer take place through convection which means no axial heat conduction. 5) Homogeneous equilibrium two-phase flow is used to describe the flow in heat exchangers with phase change. The zeotropic mixture R134a/R245fa [13] is selected as the working fluid in simulation and experiment with the following consideration: R134a and R245fa [36] are widely-used in ORC systems due to environment-friendly properties and their mixtures are attracting much attention for better system performance [9]. 2.1. Heat exchanger model based on finite volume method ORC is generally equipped with two heat exchangers, such as evaporator and condenser. Since the heat transfer process in them is similar, they can be reasonably modeled using same method. For the sake of brevity, the modeling process of evaporator is only introduced and condenser is not presented any more. The evaporator is a key system component as the place where heat absorption of fluid occurs. For simplicity and universality, this paper selects the counter-current thimble heat exchanger. The whole heat exchanger is divided into several cells along the flowing direction of working fluid, as shown in Fig. 2. The central difference scheme is implemented to connect the cell and boundary values. Specifically, the cell specific enthalpy is the arithmetic mean of specific enthalpy at the inlet and outlet. Both mass and energy conservation need to be considered in cells for working fluid and heat source. The energy conservation of the tube wall is also supposed to be considered due to its heat storage effect. The mass conservation equation of cell is:

m_ in  m_ out ¼

dm dr ¼V : dt dt

(1)

m_ in and m_ out are the inlet and outlet mass flowrate for each cell, mis the reserved mass, V represents the control volume of cell, r is the average fluid density of cell. Based on thermodynamic principles, the pressure and the specific enthalpy can be used to uniquely determine the fluid state whether in the single phase region or the two phase region.

3

Therefore, these two parameters are selected as state variables, like other authors [22,23]. The average density is approximately equal to the function of average specific enthalpy and pressure (rzrðp; hÞ). The partial derivatives of density vr=vpjh and vr=vhjp can be obtained from REFPROP v9.0.

dr vr dp vr dh þ j : ¼ j dt vp h dt vh p dt

(2)

For working fluid,

m_ f ;i  m_ f ;iþ1 ¼ Vf ;i

! vrf ;i dpf vrf ;i dhf ;i þ : jh jp vpf dt vhf ;i dt

(3)

For heat source,

m_ hs;iþ1  m_ hs;i ¼ Vhs;i

vrhs;i vhhs;i

jp

dhhs;i : dt

(4)

Subscript f and hs represent working fluid and heat source, respectively. The subscript i denotes the cell number, e.g. 1, 2, …, N. Since the pressure of the heat source is assumed to be constant, the term dphs =dt can be neglected. The energy conservation equation of cell is

dp dðmhÞ ¼ m_ in hin  m_ out hout þ Q_ þ V dt dt

(5)

hin and hout represent the specific enthalpy at the inlet and outlet. Q_ is convective heat transfer, and we specify that it is positive for absorption and negative for release. For working fluid,

  dpf m_ f ;i hf ;i  m_ f ;iþ1 hf ;iþ1 þ af ;i Afw;i Tw;i  T f ;i þ Vf ;i dt   drf ;i dhf ;i þ rf ;i : ¼ Vf ;i hf ;i dt dt

(6)

For heat source,

  m_ hs;iþ1 hhs;iþ1  m_ hs;i hhs;i þ ahs; i Ahsw;i Tw; i  T hs;i   dhhs; i drhs;i þ rhs;i : ¼ Vhs;i hhs;i dt dt

(7)

af and ahs are the corresponding convection heat transfer coefficients, which represent the value of working fluid to tube wall and heat source to tube wall, respectively. Their calculation is conducted based on correlations provided in Table 1. Since both heat source and cooling water is always in single phase, they can utilize same correlation for different cells. Respecting that working fluid undergoes phase change in heat transfer process it is

Fig. 2. Scheme of evaporator based on finite volume method.

4

J. Cai et al. / Energy 197 (2020) 117003

Table 1 Corresponding heat transfer correlations. Component

Fluid

State of phase

Correlations

Evaporator

Heat source

Single phase [37]

Nuf ¼ 0:027Re0:8 f Prf

Working fluid Working fluid

Single phase [38] Two phase [39]

1=3

ðhf =hw Þ0:14

atp ¼

Nuf ¼

ðf =8ÞðRef  1000ÞPr

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðF al Þ2 þ ðSapool Þ2

 0:35  r F ¼ 1 þ Prl l  1

rv

1 þ 12:7ðf =8Þ1=2 ððPr2=3 Þ  1Þ S ¼ ð1 þ 0:055F 0:1 Re0:16 Þ1 1 0:4 al ¼ 0:023Re0:8 l Prl ðll =dÞ

apool ¼ 55p0:12 q2b ðlgpr Þ0:55 M r Condenser

Working fluid

Single phase

Working fluid

Two phase [40]

Nuf ¼

ðf =8ÞðRef  1000ÞPr 1 þ 12:7ðf =8Þ1=2 ððPr2=3 Þ  1Þ

atp ¼ al ð1  xÞ0:8 þ

3:8x0:76 ð1  xÞ0:04

!

p0:38 r

0:4 a1 ¼ 0:023Re0:8 l Prl ðl1 =dÞ

Cooling water

Single phase

1=3

Nuf ¼ 0:027Re0:8 f Prf

reasonable to determine the phase state firstly and select the suitable correlation accordingly. Af ;w and Ahs;w are the relevant heat transfer area on the tube wall. The energy conservation equation for wall is:

    dT rw Vw cpw w;i ¼ af ; i Afw;i T f ;  Tw;i þ ahs; i Ahsw;i T hs;i  Tw; i i dt (8) There rw , Vw , cpw are density, volume and specific heat capacity of the tube wall, respectively. Therefore, five equations need to be solved to obtain the state variation in each cell. The equations from N cells can be reformulated to a set of ordinary differential equations with dimensions of 5 N. It can be calculated with the help of the solver built in SIMULINK. 2.2. Heat exchanger model based on moving boundary method The evaporator normally can be divided into three regions according to the phase state: subcooled region, two-phase region and superheated region, as shown in Fig. 3. This phenomenon gives a hint to the proposal of moving boundary method, which aims to dynamically track the lengths of the different regions. Herein, the MB model is only analyzed when three phase regions present with ignoring the condition of phase vanishing and forming. Similarly, differential equations in relation to mass and energy conservation are adopted to catch the state variation in each region and integrated along the flow direction. Since the specific derivation

ðhf =hw Þ0:14

process is shown in former literature [24,41], it is not discussed in detail in this paper. The mass conservation differential equation is:

vAr vm_ ¼ 0: þ vt vt

(9)

The energy conservation differential equation is:

_ vðArh  ApÞ vmh þ ¼ aAw ðTw  TÞ: vt vt

(10)

The energy conservation equation of the tube wall is the same as Eq. (8). The average specific enthalpy and pressure are used to determine the state of the fluid in three regions, other parameters such as average density and average temperature are their functions. In single phase region such as the subcooled and superheated, the average specific enthalpy is the arithmetic mean of the specific enthalpy at inlet and outlet. While in the two-phase region, mean void fraction is introduced to determine the state.

r ¼ rl ð1  gÞ þ rv g: h¼

rl ð1  gÞhl þ rv ghv : r

(11)

(12)

Therein, g is the void fraction which is defined as the ratio of the vapor area to the entire cross-sectional area, g is the integral average value along the direction of the tube length. Subscript l and

Fig. 3. Scheme of evaporator based on moving boundary method.

J. Cai et al. / Energy 197 (2020) 117003

v represent the saturated liquid and vapor phase, respectively. The equation recommended by Jensen [24] is adopted in the absence of precise void fraction for zeotropic mixture,

    1 þ ðrv =rl Þ2=3 2=3ln rv =rl  1 g¼1  : 2 ððrv =rl Þ2=3  1

(13)

Unlike pure fluids, the temperature of zeotropic mixtures in two-phase region varies with vapor quality. The lumped temperature should be calculated by pressure and average specific enthalpy rather than regarded as the mean of liquid temperature and vapor temperature. 2.3. Pump and turbine models Steady-state description is good choice to simplify these two components. The commonly used hydraulic diaphragm metering pump is selected as working fluid pump in this study. An algebraic expression is used to describe the pump model [24].

m_ p ¼ hv rp Vp u:

(14)

Here, hv is the volumetric efficiency, which is set to a constant of 0.9 in this paper. rp is density of the working fluid at inlet of pump, Vp is the displacement volume per revolution and u is rotation speed of the drive motor. Working fluid actually undergoes non-isentropic compression in the pump due to friction and leakage, the isentropic efficiency is set to 0.8 for simplicity. The outlet specific enthalpy can be obtained as follows:

hp;out ¼ hp;in þ

hps;out  hp;in

hp

:

(15)

hp;out and hps;out are actual non-isentropic and ideal isentropic outlet specific enthalpy, hp;in is inlet specific enthalpy and hp is isentropic efficiency. The turbine is also a key component of ORC for deciding the output power. This paper doesn’t intend to analyze the complex flow process in the turbine, and hence the flow characteristics of turbine is simply described by an algebraic equation. Generally, the actual pressure ratio between inlet and outlet of turbine is greater than critical pressure ratio resulting the supersonic flow [41e43]. Its mass flowrate can be calculated through below equation.

pffiffiffiffiffiffiffiffiffi m_ t ¼ Cv rt pt :

(16)

Cv is flow coefficient and is assumed to be a constant obtained from the design condition. The isentropic efficiency of turbine is assumed to be a constant of 0.7.

  ht;out ¼ ht;in  ht;in  hts;out ht :

(17)

ht;out and hts;out are actual non-isentropic and ideal isentropic outlet specific enthalpy, ht;in is the inlet specific enthalpy, ht is the isentropic efficiency. 2.4. Liquid receiver model The liquid receiver is placed between condenser and pump in order to reduce the fluctuation in condensation pressure and maintain the liquid state. In addition, the liquid receiver is used to charge fluid into system and its volume is generally large enough to store liquid. The dynamic model of the liquid receiver is established

5

according to the work of Eldredge [44]. The vapor and liquid in receiver are regarded as saturated state according to thermodynamic equilibrium assumption. The receiver is simplified as a simple cylindrical control volume which is insulated from the surroundings. The supply fluid coming from condenser can be either subcooled, saturated or two-phase while the exhaust fluid is always regarded as saturated liquid. The initial conservation of mass and energy (Eqs. (1) and (5)) are adopted to model receiver. The derivation results are based on two state variables: liquid level (Ll ) and receiver pressure (pr ).

  dLl dr dr dpr þ Ab Ll l þ Lv v : dt dp dp dt

(18)

dL m_ in hin  m_ out hout ¼ Ab ðrl hl  rv hv Þ l þ dt 

 dh dr dhv dr dpr þ hv Lv v  Lr : Ab rl Ll l þ hl Ll l þ rv Lv dp dp dp dp dt

(19)

m_ in  m_ out ¼ Ab ðrl  rv Þ

Lr is defined as the total height of receiver, Lv is the vapor level, Ab is the bottom area of receiver. 2.5. Simulation problem and solution 2.5.1. Flow reversal between condenser and receiver The ORC model is formed by combining the component models mentioned above. Originally, the liquid receiver pressure is always regarded as equal to the condensation pressure [44,45] whether in FV-M or MB-M. In other words, the receiver model is coupled into the condenser model and hence, the mass flowrate between them is treated as an internal state variable. However, the simulation of system model always fails when the mass flowrate of working fluid changes dramatically. Among those failure cases, the fluid mass flowrate between condenser and receiver is usually found to be negative. The undesirable flow reversal is only the result of mathematical calculation and doesn’t occur in the actual practice. This abnormal situation could cause collapse of the condenser model based on central difference scheme and result in simulation failures. The flow reversal could be explained qualitatively. When the pressure in receiver is always kept equal to that in condenser, the sudden variation of pressure in condenser will cause dramatic change of fluid enthalpy in receiver due to the saturation assumption. The update of fluid enthalpy in receiver only relies on convective energy flow. However, the enthalpy of fluid coming from condenser changes more slowly than pressure in receiver. Thus the inlet flowrate will turn negative to realize the rapid update of fluid energy in receiver. Similar generalized demonstration of the flow reversal can be found in the work of Schulze [46] and Quoilin [47]. In order to maintain a positive flow between these two components in simulation, a valve model obtained from Ref. [45] is introduced to keep some pressure difference between them.

m_ v ¼ Xv Av

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð2rv Dpv Þ:

(20)

Av is the cross-sectional area of valve when fully open, Xv is the actual opening degree. rv is the inlet density and Dpv is pressure difference between inlet and outlet. Under the effect of the pressure difference, the fluid will always flow from condenser to receiver. 2.5.2. Chattering in FV-M This paper focuses on the subcritical ORC, which is characterized by phase transition occurring in both evaporator and condenser. The phase transition brings deterioration to the robustness and

6

J. Cai et al. / Energy 197 (2020) 117003

derivative smoothing method as recommended in Ref. [47] to circumvent the problem and promote simulation robustness. The idea is to smooth out the density derivative discontinuity using a linear interpolation function in the transition region which is distinguished by vapor quality x (x2½0; d, d is the user-defined scale of vapor quality for interpolation). Since the heat transfer coefficients of working fluid in three phase regions are calculated by different correlations, the discontinuity between them should not be ignored. The discontinuity could result in abnormal energy flow in convection heat transfer between fluid and tube when state switch occurs in cell. A trigonometric function [45] is introduced to conduct the interpolation in the region near the saturated liquid point and the saturated vapor point.

 Fig. 4. Original chattering and linearized improvement in FV-M.

atp;x ¼

0:5ðal þ aε Þþ0:5ðal  aε Þcosðpx=εÞfor0
accuracy of the dynamic model [28]. Chattering could be often observed in the simulation based on FV-M, as illustrated in Fig. 4. The original chattering situation of MORC using R134a/R245fa (0.5/ 0.5) based on FV-M is demonstrated in Fig. 4, it is obtained under the step change of working fluid mass flowrate. This phenomenon is attributed to the discontinuity of variables in simulation process. The partial derivatives of density with pressure vr= vpjp and with

enthalpy vr=vhjh varies with vapor quality, the situation of R134a/ R245fa (0.5/0.5) is displayed as an example in Fig. 5. The original data stands for those directly calculated based on REFPROP v9.0 and the maximum discontinuity can be found near the saturated liquid point. When thermodynamic state of fluid in cell varies, the discontinuity of two partial derivatives gives rise to a spurious numerical internal mass flow in control volume. This can greatly reduce the simulation speed and increase numerical errors, eventually result in simulation failure. This paper chooses the density

Therein, atp;x is convection heat transfer coefficient of fluid at vapor quality of x. ε is the user-defined scale of vapor quality for interpolation. In the present work, both d and ε are imposed to 0.1. With the recommended linearization, the mutations of parameters will be eliminated. And the possibility of the calculated variables exceeding the acceptable boundaries is greatly reduced. Although these measures will slightly deteriorate the accuracy of the model, the author appreciates the positive effect on avoiding chattering and speeding up the simulation. The improvement based on linearized method is represented by dot line in Fig. 4. It can be seen that the magnitude of the abnormal chatter is greatly reduced.

3. Test bench description A kW-scale ORC test bench is established in order to verify the

Fig. 5. Original and linearized partial derivative of density.

! vr vr vr vr x ¼ þ  : vp h;x vp h;0 vp h;d vp h;0 d

(21)

  vr vr vr vr x ¼ þ  : vh p;x vh p;0 vh p;d vh p;0 d

(22)

J. Cai et al. / Energy 197 (2020) 117003

7

Fig. 6. Schematic layout of experiment system.

Fig. 7. Photograph of ORC test bench.

accuracy of proposed model. The test bench is consisted of three parts which are heat source, working fluid and heat sink, respectively. Schematic diagram of the experiment system is given in Fig. 6 and the photograph of real test bench is given in Fig. 7. Heated air is used as heat source of ORC system because of low cost and easy access. The equipped air heater can produce adjustable high temperature air. The flowrate and temperature of the heated air can be controlled by the operating frequency of the fan and the power of the heating resistance wire. The cooling water is selected as the heat sink. The working fluid undergoes a series of typical flow and physical change. Firstly, liquid fluid coming from the liquid receiver is pumped to the evaporator and vaporized by high temperature air. Then, the vapor from evaporation passes through the expansion valve and subsequently is cooled into liquid in the condenser.

Finally, the fluid returns to the receiver. The receiver is selfdesigned cylindrical container with the volume of 5 L and a hydraulic diaphragm metering pump is purchased from Milton Roy. The evaporator is designed and processed independently, whose inner pipe is spiraled into air conduit aiming to get larger heat transfer area. Working fluid flows through the inner tube and air flows through the outer tube. Since the turbine of ORC is under delivering, a self-made expansion valve is employed to temporarily take the place of the turbine. The condenser is a brazed-plate heat exchanger purchased from SWEP for high integration and efficiency. The damper and the steady-flow compensator are arranged between the pump and the flowrate meter in order to reduce the flow fluctuation and improve the measurement accuracy of flowrate. The test bench has also been equipped with measuring

8

J. Cai et al. / Energy 197 (2020) 117003

Table 2 Measuring instruments and accuracy of main parameters. Parameter

Sensor

Scale

Accuracy

Temperature Pressure

Sheathed thermocouple sensors with first-class precision Differential pressure transmitter (air) Differential pressure transmitter (fluid) Air flowmeter Coriolis mass flowmeter Electromagnetic flowmeter

60e650
C 0e0.5 MPa 0e4 MPa 0e500 Nm^3/h 0e100 kg/h 0e5 m^3/h

±0.4% ±0.065% ±0.065% ±1.5% ±0.2% ±0.5%

Flowrate

instruments, such as pressure sensors, temperature sensors and flowmeters. Parameters about type, scale and accuracy of corresponding instruments are listed in Table 2. Data acquisition of measurement is by means of the multi-channel input modules and the signal is transmitted to the computer. This data acquisition system can achieve the recording speed of 1 Hz. The tubes and components of the whole test bench are all wrapped in thermalinsulation materials to minimize the impact of environmental conditions. The uncertainty calculation of indirect measurement parameters is conducted based on the formula recommended by Ref. [48].

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u n  2 uX vy dy ¼ t dxi 2 : vx i¼1

(24)

where parameter y is the function of the measurements (x1 , x2 , …, xn ), dy and dxi are uncertainty of y and xi , respectively. For the sake of demonstration, the uncertainty of heat absorption as a key parameter of ORC system is introduced. The steady heat absorption of fluid is generally calculated by Eq. (25).

_ out ðTout ; pÞ  hin ðTin ; pÞ: Q_ ¼ m½h

(25)

And the calculation of uncertainty is conducted by Eq. (26) below.

Fig. 8. Trajectories of measured, simulated and filtered mass flowrate of R134a/R245fa (0.45/0.55).

changes take place on mass flowrate of working fluid while the input parameters of heated air and cooling water almost remained unchanged. The frequency of pump drive motor undergoes a series of step change (45-40-45-50-45 Hz), the whole experiment lasts for

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 12 3 u 12 0 20 !2 3 u  2 u vhout vhin 5 2 vh vhin C 7 2 2 2 2 dQ_ ¼ u þ dT þ m_ 2 4@ out A þ B @ A 5dp : tðhout  hin Þ dm_ þ m_ vTout vTin vp vp

Based on steady state experiment data, the maximum relative uncertainty concerning heat absorption of working fluid is 1.18%. This indicates the test bench is of high reliability and accuracy. 4. Model validation Two kinds of dynamic model of MORC have been established based on above work. In this section, both FV-M and MB-M are validated by the experimental data from the test bench since there are no verification research available on MORC dynamic model in the literature. In consideration of trade-off between simulation accuracy and speed, the number of exchanger cells is selected to be 20 for FV-M, referring to other authors [22,35]. In order to start MBM effectively, three phase regions are present in evaporator and condenser at the beginning of the test. The working fluid charged in test bench is R134a/R245fa (0.45/0.55) which is obtained by weighing method. During the whole experiment, a series of step

(26)

a total of 2760 s and transition points are 343, 932, 1558 and 2132 (s), respectively. It should be pointed out that both FV-M and MB-M are corrected and initialized with experimental data at the initial state. The main correction method is to keep the length of the heat exchange tube same as that of the test bench by multiplying the original heat exchange coefficients by some factors. Trajectories of measured, simulated and filtered mass flowrate of R134a/R245fa (0.45/0.55) is shown in Fig. 8. The simulated is calculated from pump model and the filtered is obtained from the measured through data processing. Obviously, the simulated is in poor agreement with the measured, which is attributed to the treatment of constant volumetric efficiency of pump in model. The actual volumetric efficiency of pump varies with working condition. However, the improvement is difficult by modifying the volumetric efficiency in the absence of pump performance curve. Thus, adopting filtered mass flowrate is a good option for maintaining the consistency in boundary condition of simulation and

J. Cai et al. / Energy 197 (2020) 117003

9

Fig. 9. Trajectories of some MORC system parameters obtained from experiment and simulation.

Table 3 Calculated relative differences of some parameters comparing to experiment. FV-M

Evaporation pressure Condensation pressure Outlet temperature of fluid Outlet temperature of air

MB-M

Dxmax;rel

Dxrel

Dxmax;rel

Dxrel

4.92% 1.73% 10.95% 1.07%

1.54% 0.66% 3.60% 0.41%

4.84% 1.76% 13.95% 1.30%

1.31% 0.54% 4.13% 0.46%

experiment. Herein, the maximum and average relative errors of the two models with respect to experiment are calculated by Eq. (27) and Eq. (28), respectively.

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! ð~xk  xk Þ2 Dxmax;rel ¼ max ; k2½0; ktot : ~x

(27)

k

pffiffiffiffiffiffiffiffiffiffiffiffiffiffi!

Dxrel ¼

tot Skk¼0

ð~xk xk Þ2 ~x k

ktot

:

(28)

xk indicates the calculated value at time point k, while ~xk is the corresponding measured value and total number of time points is ktot .

The corresponding results of experiment and simulation concerning system main parameters are shown in Fig. 9. The simulation of MB-M was interrupted in 1636 s, this resulted from the superheated region vanishing. In contrast, the FV-M was valid in the whole process. Simulation results of two models generally fitted the variation trend of experiment from the view of Fig. 9. Table 3 presents the errors about four parameters and all the mean relative errors are apparently less than 5%. This result indicates that two system models could achieve similar high simulation accuracy. In addition, the high model accuracy indicates that ignoring composition shift in dynamic modeling of MORC is acceptable. Nevertheless, there are some distinct differences between simulation and experiment. The prediction ability of two models deteriorates when the working condition deviating from the initial state, particularly reflected in evaporation pressure and outlet temperature of fluid. As shown in Table 3, the maximum relative difference about evaporation pressure is around 5%, the value of fluid outlet temperature is 11% for FV-M (14% for MB-M). This is mainly caused by uncertainties in heat transfer correlations due to complexity of flow and inaccuracy of the Nusselt correlation [19]. The correction factor based on the initial steady state actually has been applied to the heat transfer coefficient, but apparently its effect is limited. What’s more, the adiabatic assumption makes simulation in not accord with reality. Another cause is attributed to the constant discharge coefficient which is actually variable under different working condition. Fig. 9 also indicates that the simple MB-M established in this

10

J. Cai et al. / Energy 197 (2020) 117003

Table 4 Main system parameters of ORC model. Parameter

Value

Unit

Inlet temperature of heated air Mass flowrate of heated air Inlet temperature of cooling water Evaporation pressure Superheated degree Condensation temperature Subcooled degree Isentropic efficiency of turbine Isentropic efficiency of pump Inside diameter of fluid tube Outside diameter of fluid tube Inside diameter of air tube Inside diameter of water tube Total length of evaporator Total length of condenser

250 0.08 15 3 10 25 5 0.7 0.8 0.006 0.008 0.155 0.014 18 14


C kg/s
C MPa
C
C
C e e m m m m m m

paper is helpless under the abnormal condition like disappear of superheated region. There are already improved methods like switching moving boundary to deal with phase vanishing and forming, however, the derivation of switching models [25,28] is quite complex and still subject to the void coefficient which is uncertain for zeotropic mixture. In contrast, the FV-M is applicable to the situation of phase vanishing and forming. Therefore, the following analysis about MORC is based on the FV-M. 4.1. Model application-dynamic analysis for MORC This section focuses on the dynamic process of MORC system with mixture at different mass fraction. As ORC system contains many degrees of freedom, some variables should be kept same in order to stay on the same benchmark at the initial state. Those specified parameters are shown in Table 4. Response characteristics are characterized by variation direction and response speed. Here, pure R134a and R245fa are nominally named as R134a/R245fa (1.0/ 0.0) and R134a/R245fa (0.0/1.0) like their mixtures for a unified fluid description. Five proportions in mass fraction are selected for research: 0.0/1.0, 0.2/0.8, 0.5/0.5, 0.8/0.2 and 1.0/0.0. 4.2. Response behavior presentation The basic ORC for exhaust waste heat recovery is usually equipped with four major components: evaporator, condenser, pump and turbine. Among them, pump and turbine response fast for small capacity. Condenser is accompanied with small

temperature difference and high heat transfer coefficient and this results in faster response and smaller variation comparing to evaporator. Specifically, evaporator plays a key role on the system performance. Hence the dynamic behavior of evaporator is chosen to represent the system situation. In this paper, evaporation pressure and outlet specific enthalpy of working fluid in evaporator are the selected parameters to show the response situation. The following enthalpy is referred to the relative deviation from the initial value for getting same base benchmark. The simulations are arranged as follows: all dynamic simulation lasts for 500 s, they start from the design point and remain steady until 50 s at which a step disturbance is imported, then these models run freely over the remaining 450 s. Considering that exhaust frequently fluctuates and fluid mass flowrate is easily adjusted, their disturbances are selected as objects of investigation. Specifically, the working fluid flowrate (mf) decreases by 10%, the exhaust mass flowrate (mg) increases by 10%, and the inlet exhaust temperature (Tg) decreases by 20 K. Figs. 10, 11 and 12 describe the response results under those three disturbances, respectively. When the working fluid flowrate decreases, the pressure in the evaporator decreases and the outlet specific enthalpy of fluid increases. It can be explained by that the outlet mass flowrate of fluid is higher than the inlet for a short period after the step due to flow inertia and this causes the accumulative quantity of working fluid in evaporator decreases. The balance between convective heat flow and heat transfer energy flow is reached at the initial steady state, the sharp decrease of fluid mass flowrate leads to decrease in convective and heat transfer energy flow, and the decline of convective energy flow is more. Therefore, the specific enthalpy of fluid at the outlet increases and its temperature also goes up. The disturbance in temperature and mass flowrate of exhaust influences the evaporation of working fluid mainly through heat transfer performance. The increase of exhaust flowrate enhances heat transfer and promotes evaporation thus leading to higher evaporation pressure, while the decrease of inlet exhaust temperature weakens heat transfer and reduces evaporation thus resulting in lower evaporation pressure. The variation direction of outlet specific enthalpy of working fluid can also be explained by the balance between convective energy flow and heat transfer energy flow. These figures all reveal that the variation trend of zeotropic mixtures at different mass fraction is same under the identical input of disturbance, the variation amplitude and response speed are different though.

Fig. 10. Dynamic response behavior when working fluid mass flowrate decreases.

J. Cai et al. / Energy 197 (2020) 117003

11

Fig. 11. Dynamic response behavior when exhaust mass flowrate increases.

Fig. 12. Dynamic response behavior when exhaust temperature decreases.

Fig. 13. Parameters of response characteristics under three kinds of disturbances.

4.3. Response speed comparison Two performance indicators calculated from dynamic process are adopted to evaluate the response speed of the system: rise time and time constant [41]. Rise time refers to the required time to go from 10% to 90% of the final value. Time constant indicates the time required to reach 1  1=ez63:2% of the final value. Typically, the smaller these two values are, the faster the system responds. Fig. 13

shows the calculated results of zeotropic mixtures at different mass fraction. It can be seen that with these three disturbances, time indicators of pure R134a is always smaller than pure R245fa. This result is consistent with Shu’s [41] research on pure fluid. As for mixtures, rise time and time constant generally have a tendency to decrease with the increase of mass fraction of R134a which represents the component having fast respond ability. It can be explained by that physical properties of mixtures move towards the

12

J. Cai et al. / Energy 197 (2020) 117003

increasing pure component. It can be seen that response speed of zeotropic mixture is generally between two pure components and so that response speed of mixtures can be roughly predicted according to this principle. Obviously, response speed of evaporation pressure is faster than outlet specific enthalpy and it can be concluded that the slow response speed of the system is mainly determined by the enthalpy of working fluid. The rough physical interpretation [43] is that the faster hydraulic disturbance propagates at the speed of sound and the slower thermal disturbance propagates at the speed of flow. For evaporation pressure, response speed under fluid flowrate disturbance is obviously faster than other two disturbances. For fluid outlet enthalpy, response speed under three kinds of disturbances is similar.

5. Conclusion Dynamic models of ORC using zeotropic mixture are established by finite volume method and moving boundary method, respectively. These two models are well verified by experiment. The dynamic response characteristics of MORC under three kinds of disturbances are studied based on FV-M. Corresponding results can be summarized as follows: (1) While coupling the condenser and receiver model, the flow reversal between them could be encountered. A feasible method to solve the problem is to add a valve model between them. Moreover, the linearization of partial derivative of density and heat transfer coefficient is recommended to eliminate the chattering in simulation of FV-M. (2) FV-M and MB-M have similar high simulation accuracy for MORC according to the experimental validation. Therefore, it proves that composition shift in MORC could be neglected and the average void fraction expression of pure fluid also shows effectiveness to mixture in MB-M. These efforts could provide confidence for the following studies on the dynamic progress of MORC. (3) The zeotropic mixtures with different mass fraction have the similar variation trend under same disturbances. The response speed of system with R134a/R245fa increases with the mass fraction of R134a. The response speed of evaporation pressure is faster than outlet specific enthalpy. Specifically, the rise time and time constant of evaporation pressure under fluid mass flowrate disturbance are the smallest comparing to disturbance in exhaust mass flowrate and temperature. This paper extends the dynamic ORC research from pure fluids to zeotropic mixtures. The system control for MORC will be effectively developed based on the proposed model in the future study. Actually, the experimental validation of models reveals that there are some distinct difference between simulation and experiment. We will study this in detail and try to make improvements.

Declaration of competing interest The authors declare no conflict of interest.

Acknowledgment This work was supported by a grant from the National Natural Science Foundation of China (No. 51676133).

Nomenclature

Symbols m_ m t V p h Q_

Mass flowrate (kg/s) Reserved mass (kg) Time (s) Volume (m3) Pressure (kPa). Specific enthalpy (kJ/kg).

Qe Qa Qc Qcw

Heat transfer rate (kW). Area (m2) Temperature (K) Specific heat capacity (J/kg$K) Stodola’s coefficient () Isentropic efficiency of turbine () Length (m) Opening degree of valve () Vapor quality () Time point () Mass flowrate of working fluid (kg/s) Mass flowrate of exhaust (kg/s) Inlet temperature of exhaust (K) Evaporation pressure when mf decreases by 10% Evaporation pressure when mg increases by 10% Evaporation pressure when Tg decreases by 20 K Outlet enthalpy of evaporator when mf decreases by 10% Outlet enthalpy of evaporator when mg increases by 10% Outlet enthalpy of evaporator when Tg decreases by 20 K Isentropic efficiency of pump () Variable () Heat adsorption of fluid in evaporator Heat release of air in evaporator Heat release of fluid in condenser Heat adsorption of water in condenser

Subscripts in out f hs i fw w hsw l v p t r b max rel

Inlet Outlet Working fluid Heat source Cell number Fluid to tube wall Tube wall Heat source to tube wall Liquid at saturated sate Vapor at saturated state/Valve Pump Turbine Receiver Bottom of receiver Maximum Relative

A T cp Ks

ht L Xv x k mf mg Tg Pe_mf Pe_mg Pe_Tg heout_mf heout_mg heout_Tg

hp x

Abbreviations ORC Organic Rankine Cycle FV-M Model based on Finite Volume method MB-M Model based on Moving Boundary method MORC Organic Rankine Cycle using Mixture References [1] Nguyen TQ, Slawnwhite JD, Boulama KG. Power generation from residual industrial heat. Energy Convers Manag 2010;51:2220e9. [2] Shi L, Shu G, Tian H, Deng S. A review of modified Organic Rankine cycles

J. Cai et al. / Energy 197 (2020) 117003

[3]

[4]

[5] [6] [7] [8]

[9]

[10]

[11]

[12]

[13] [14]

[15] [16]

[17] [18]

[19]

[20]

[21]

[22]

[23]

[24]

(ORCs) for internal combustion engine waste heat recovery (ICE-WHR). Renew Sustain Energy Rev 2018;92:95e110. €m H-E. A review of turbocompounding as a waste heat Aghaali H, Ångstro recovery system for internal combustion engines. Renew Sustain Energy Rev 2015;49:813e24. Delgado-Torres AM, García-Rodríguez L. Comparison of solar technologies for driving a desalination system by means of an organic Rankine cycle. Desalination 2007;216:276e91. s C, Kaparaju P. Geothermal energy: power plant technology and Moya D, Alda direct heat applications. Renew Sustain Energy Rev 2018;94:889e901. Zhang X, He M, Wang J. A new method used to evaluate organic working fluids. Energy 2014;67:363e9. Miao Z, Zhang K, Wang M, Xu J. Thermodynamic selection criteria of zeotropic mixtures for subcritical organic Rankine cycle. Energy 2019;167:484e97. Liu C, Gao T. Off-design performance analysis of basic ORC, ORC using zeotropic mixtures and composition-adjustable ORC under optimal control strategy. Energy 2019;171:95e108. Sadeghi S, Maghsoudi P, Shabani B, Gorgani HH, Shabani N. Performance analysis and multi-objective optimization of an organic Rankine cycle with binary zeotropic working fluid employing modified artificial bee colony algorithm. J Therm Anal Calorim 2019;136:1645e65. Bamorovat Abadi G, Kim KC. Investigation of organic Rankine cycles with zeotropic mixtures as a working fluid: advantages and issues. Renew Sustain Energy Rev 2017;73:1000e13. Tian H, Wu M, Shu G, Liu Y, Wang X. Experimental and theoretical study of flammability limits of hydrocarboneCO2 mixture. Int J Hydrogen Energy 2017;42:29597e605. El-Sayed AR, El Morsi M, Mahmoud NA. A review of the potential replacements of HCFC/HFCs using environment-friendly refrigerants. Int J AirCond Refrig 2018;26. Abadi GB, Yun E, Kim KC. Experimental study of a 1 kw organic Rankine cycle with a zeotropic mixture of R245fa/R134a. Energy 2015;93:2363e73. Shu G, Yu Z, Tian H, Liu P, Xu Z. Potential of the transcritical Rankine cycle using CO2-based binary zeotropic mixtures for engine’s waste heat recovery. Energy Convers Manag 2018;174:668e85. Wang JL, Zhao L, Wang XD. A comparative study of pure and zeotropic mixtures in low-temperature solar Rankine cycle. Appl Energy 2010;87:3366e73. Shu G, Li X, Tian H, Shi L, Wang X, Yu G. Design condition and operating strategy analysis of CO 2 transcritical waste heat recovery system for engine with variable operating conditions. Energy Convers Manag 2017;142:188e99. Liu P. Engine load effects on the energy and exergy performance of a medium cycle/organic rankine cycle for exhaust waste heat recovery. Entropy 2018;20. Benato A, Kærn MR, Pierobon L, Stoppato A, Haglind F. Analysis of hot spots in boilers of organic Rankine cycle units during transient operation. Appl Energy 2015;151:119e31. Horst TA, Rottengruber H-S, Seifert M, Ringler J. Dynamic heat exchanger model for performance prediction and control system design of automotive waste heat recovery systems. Appl Energy 2013;105:293e303. Wang X, Shu G, Tian H, Liu P, Li X, Jing D. Engine working condition effects on the dynamic response of organic Rankine cycle as exhaust waste heat recovery system. Appl Therm Eng 2017;123:670e81. Quoilin S, Aumann R, Grill A, Schuster A, Lemort V, Spliethoff H. Dynamic modeling and optimal control strategy of waste heat recovery Organic Rankine Cycles. Appl Energy 2011;88:2183e90. Ni J, Zhao L, Zhang Z, Zhang Y, Zhang J, Deng S, Ma M. Dynamic performance investigation of organic Rankine cycle driven by solar energy under cloudy condition. Energy 2018;147:122e41. Yousefzadeh M, Uzgoren E. Mass-conserving dynamic organic Rankine cycle model to investigate the link between mass distribution and system state. Energy 2015;93:1128e39. Jensen JM. Moving boundary models for dynamic simulation of two-phase

13

flows. In: The second international modelica conference, Germany; 2002. [25] Huster WR, Vaupel Y, Mhamdi A, Mitsos A. Validated dynamic model of an organic Rankine cycle (ORC) for waste heat recovery in a diesel truck. Energy 2018;151:647e61. [26] Wang X, Shu G, Tian H, Liu P, Jing D, Li X. Dynamic analysis of the dual-loop Organic Rankine Cycle for waste heat recovery of a natural gas engine. Energy Convers Manag 2017;148:724e36. [27] Wei D, Lu X, Lu Z, Gu J. Dynamic modeling and simulation of an Organic Rankine Cycle (ORC) system for waste heat recovery. Appl Therm Eng 2008;28:1216e24. [28] Vaupel Y, Huster WR, Holtorf F, Mhamdi A, Mitsos A. Analysis and improvement of dynamic heat exchanger models for nominal and start-up operation. Energy 2019;169:1191e201. [29] Zhao L, Bao J. The influence of composition shift on organic Rankine cycle (ORC) with zeotropic mixtures. Energy Convers Manag 2014;83:203e11. [30] Zhou Y, Zhang F, Yu L. The discussion of composition shift in organic Rankine cycle using zeotropic mixtures. Energy Convers Manag 2017;140:324e33. [31] Zhang J, Mondejar ME, Haglind F. General heat transfer correlations for flow boiling of zeotropic mixtures in horizontal plain tubes. Appl Therm Eng 2019;150:824e39. [32] Xu B, Rathod D, Yebi A, Filipi Z, Onori S, Hoffman M. A comprehensive review of organic rankine cycle waste heat recovery systems in heavy-duty diesel engine applications. Renew Sustain Energy Rev 2019;107:145e70. [33] MATLAB R2015b. https://ww2.mathworks.cn. [34] REFPROP v9.0. https://www.nist.gov/srd/refprop. [35] Chowdhury JI, Nguyen K, Thornhill D. Dynamic model of supercritical organic Rankine Cycle waste heat recovery system for internal combustion engine. Int J Automot Technol 2017;18(4):589e601. [36] Maraver D, Royo J, Lemort V, Quoilin S. Systematic optimization of subcritical and transcritical organic Rankine cycles (ORCs) constrained by technical parameters in multiple applications. Appl Energy 2014;117:11e29. [37] Sieder EN, Tate GE. Heat Trans Press Drop of Liq Tubes 1936;28:1429e35. [38] Gnielinski V. New equations for heat and mass transfer in turbulent pipe and channel flow. Int Chem Eng 1976;16:359e68. [39] Liu Z, Winterton RH. Transfer, a general correlation for saturated and subcooled flow boiling in tubes and annuli. Base Nucleate Pool Boiling Equ 1991;34:2759e66. [40] Shah MM. Transfer, a general correlation for heat transfer during film condensation inside pipes. Int J Heat and Mass Trans 1979;22:547e56. [41] Shu G, Wang X, Tian H, Liu P, Jing D, Li X. Scan of working fluids based on dynamic response characters for Organic Rankine Cycle using for engine waste heat recovery. Energy 2017;133:609e20. [42] Xu B, Rathod D, Kulkarni S, Yebi A, Filipi Z, Onori S, Hoffman M. Transient dynamic modeling and validation of an organic Rankine cycle waste heat recovery system for heavy duty diesel engine applications. Appl Energy 2017;205:260e79. [43] Jensen HTJM. Moving boundary models for dynamic simulation of two-phase flows. 2002. [44] Eldredge BD. Improving the accuracy and scope of control-oriented vapor compression cycle system models. 2006. [45] Quoilin S. Sustainable energy conversion through the use of Organic Rankine Cycles for waste heat recovery and solar applications. Belgium: University of ge; 2011. Lie €ber M, Tegethoff W. A limiter for preventing singularity in [46] Schulze C, Gra simplified finite volume methods45; 2012. p. 1095e100. [47] Quoilin S, Bell I, Desideri A, Dewallef P, Lemort V. Methods to increase the robustness of finite-volume flow models in thermodynamic systems. Energies 2014;7:1621e40. [48] Shu G, Zhao M, Tian H, Wei H, Liang X, Huo Y, Zhu W. Experimental investigation on thermal OS/ORC (Oil Storage/Organic Rankine Cycle) system for waste heat recovery from diesel engine. Energy 2016;107:693e706.