- Email: [email protected]
Dynamic behavior of supercritical organic Rankine cycle using zeotropic mixture working fluids
Journal Pre-proof Dynamic behavior of supercritical organic Rankine cycle using zeotropic mixture working fluids
Xiaoxue Chen, Chao Liu, Qibin Li, Xurong Wang, Shukun Wang PII:
S0360-5442(19)32271-6
DOI:
https://doi.org/10.1016/j.energy.2019.116576
Reference:
EGY 116576
To appear in:
Energy
Received Date:
06 August 2019
Accepted Date:
17 November 2019
Please cite this article as: Xiaoxue Chen, Chao Liu, Qibin Li, Xurong Wang, Shukun Wang, Dynamic behavior of supercritical organic Rankine cycle using zeotropic mixture working fluids, Energy (2019), https://doi.org/10.1016/j.energy.2019.116576
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.
Journal Pre-proof
Dynamic behavior of supercritical organic Rankine cycle using zeotropic mixture working fluids
Xiaoxue Chen, Chao Liu*, Qibin Li*, Xurong Wang, Shukun Wang
Key laboratory of low-grade Energy Utilization Technologies and Systems, Ministry of Education, School of Energy and Power Engineering, Chongqing University, Chongqing 400030, China *Corresponding author: 1. Chao Liu. E-mail address: [email protected] 2. Qibin Li. E-mail address: [email protected]
Abstract As a prominent technology for recovering low-grade waste heat, supercritical organic Rankine cycle (ORC) exhibits a better performance due to the higher endothermic temperature and better thermal match with the heat source. The irreversibility in system can be reduced greatly when the zeotropic mixtures are used as working fluid to match the heat source and sink profiles. Affected by the fluctuations in waste heat sources, there is a challenge for ORC to recovery waste heat. An improved dynamic model of supercritical ORC using zeotropic mixture R134a/R32 as working fluid is developed and dynamic behaviors of supercritical ORC are analyzed. It is found that an abnormal fluctuation may occur in some parameters due to the effects around pseudo-critical point and thermal inertia of the
Journal Pre-proof heater. A fitting correlation to predict the response time based on different heat transfer coefficients of heater is found. Besides, three dynamic regimes are defined to investigate the effects of heat source frequencies and system thermal inertia on the dynamic response. As the heat source frequency increases or the heat exchange in the heater is enhanced, the fluctuation amplitude of the pressure decreases in the heater. Key words: Supercritical organic Rankine cycle, zeotropic mixture, dynamic characteristic, thermal inertia, dynamic regime
1. Introduction In recent years, the issue of energy shortage and carbon dioxide emission has been increasingly serious. The utilization of renewable energy source has been attracted more attention to hedge against rising costs of fossil fuels and reduction in greenhouse gas emissions. Furthermore, there is a high potential for waste heat recovery in many industries such as cement industries, chemical industries, iron and steel industries. The organic Rankine cycle (ORC) is a superior technology for utilization of low-grade energy such as solar energy, geothermal energy and waste heat, etc. Generally, two types of ORC for different operation conditions are mainly included: subcritical and supercritical ORC. Li et al. [1] also proposed a novel ORC that couples supercritical and subcritical heat absorption processes, which can better enhance the system efficiency. The comparisons between the performance of
Journal Pre-proof supercritical and subcritical ORC have been conducted in some researches [2, 3]. Supercritical approach, by contrast, can achieve a better performance of system due to the higher endothermic temperature and better thermal match with the heat source [4]. Thermal and exergy efficiencies between the subcritical and supercritical ORC for biogas fueled combined heat and power engine system were compared by Yagli et al. [3]. Their results showed that the supercritical ORC can make much network output and a higher thermal and exergy efficiency of system. The selection of working fluid is crucial to the ORC system performance. Carbon dioxide (CO2) with relatively low critical temperatures and pressures is widely used in supercritical ORC. Some organic fluids such as hydrocarbons [5] and refrigerants R123 [6], R134a [7], R245fa, R143a [8] has also been studied as working fluid of supercritical ORC. Since the phase change of pure working fluids has a large temperature difference to result in a bad thermal match with the heat source or heat sink, the zeotropic mixtures emerge great advantages due to its temperature glide. Chen et al. [9] introduced the idea of using zeotropic mixtures as working fluid for supercritical ORC. The comparisons between the conventional ORC and supercritical ORC indicate that the higher thermal efficiency and exergy efficiency of supercritical ORC can be achieved by using a zeotropic mixture R134a/R32. Radulovic et al. [10] made a comparison on the operation performance of supercritical ORC using six zeotropic mixtures and found the thermal efficiency of system using pure working fluid is overall lower than using zeotropic mixtures. Dong et al.[11] studied zeotropic
Journal Pre-proof mixtures of siloxanes as working fluids in high-temperature ORC. The zeotropic mixture working fluids have been indicated to has a great potential in economic profitability [12]. In order to achieve a better performance, many studies have been carried out, especially for the steady state. In fact, affected by some uncertain factors, the parameters of heat source often fluctuate over the time in actual operating conditions. In addition, the operating mode switching may result in an inefficiency, particularly for the exhaust heat recovery of internal combustion engines [13]. It is significant to predict the dynamic behavior of system under the unsteady state. Quoilin et al. [14] proposed a dynamic model of the ORC focusing specifically on the time-varying performance of heat exchangers. Three different control strategies also were tested and compared to optimize the working conditions. Horst et al. [15] presented an approach for predicting dynamic Rankine Cycle performance and fuel saving potential in passenger car. A multivariable control strategy was proposed in order to achieve both transient performance and steady-state energy saving [16]. Wang et al. [17] established a dynamic model to study the dynamic behaviors and part-load performance of the double loop ORC. The effects of the mass flow rate of high temperature cooling water on the system performance are researched as well. Chowdhury et al. [18] presented a Finite Volume (FV) dynamic model of supercritical ORC and analyzed the dynamic behaviors and transient responses of the system based on the transient heat source. Due to the limitation of time-intensive calculation of FV
Journal Pre-proof model, they also developed a fuzzy-based dynamic model for real-time applications [19]. In view of seasonal and meteorological factors, Wang et al. [20] established a series of supercritical ORC thermodynamic models under the design and off-design conditions. Olumayegun et al. [21] proposed a dynamic model of supercritical CO2 power cycle for cement industry waste heat recovery and found it is preferable to allow the temperature of the turbine inlet to vary with the transient heat source. However, most present dynamic studies are focused on the subcritical system and the literatures on the dynamic analysis of supercritical ORC are seldom. In addition, the thermo-physical properties of the fluids in supercritical pressures are highly variable with the temperature and existing dynamic models cannot eliminate the effects of such variation. This leads to a loss of accuracy or computational efficiency of models. Therefore, an improved model to consider the effect of variable thermo-physical property in supercritical fluid is developed in this work. Both the model accuracy and computational efficiency can be satisfied simultaneously. The dynamic behavior of supercritical ORC is analyzed in this work. The effects around pseudo-critical point and thermal inertia of heater are investigated. The model is discretized by moving boundary, thus variable heat transfer area of a discretized cell can be used as a detection of heat transfer capacity, and furthermore the position of pseudo-critical point in the heater can be found. In addition, the dynamic responses of the supercritical system with different thermal inertia and variable frequencies of heat source are presented. A mathematical correlation between the response time and
Journal Pre-proof transfer coefficients of the heater is found, which can give guidance for control and design of the system on the basis of different demands and application conditions.
2. System description and cycle process Xu [8] pointed out that the critical temperature of most suitable working fluid is 40-65 ℃ lower than the flue gas inlet temperature and the optimum condition is supercritical. Taking the fluids’ properties and operation conditions into consideration, the mixture of R134a/R32 (mass fraction is [0.7/0.3]) is selected in this study, which has been turned out to achieve a high thermal and exergy efficiency of system [9]. The schematic layout and T-s diagram of the supercritical ORC with zeotropic mixture of 0.7R134a/0.3R32 are shown in Fig.1. The working fluid is heated by the flue gas and becomes supercritical vapor directly from the liquid in the heater (4-1). The supercritical vapor enters into the turbine and completes the expansion to generate power. Then the superheated vapor from the turbine flows into the condenser where it is cooled down to subcooled liquid by cooling water but non-isothermally (2-3). After that, the liquid is pressurized (3-4) to be supercritical fluid into the heater, and a new cycle begins.
Journal Pre-proof
Fig.1 (a) Schematic of ORC system and (b) T-S diagram of the process
The main parameters for the initial design point are listed in Table 1 [9, 22]. The critical pressure of 0.7R134a/0.3R32 is 5.13 MPa [9]. The high pressure of initial designed cycle is set 7 MPa and the inlet temperature of working fluid in turbine is set 413 K, which is 10 K lower than the temperature of flue gas [9]. The droplets formation can be avoided in the expansion process [23]. The condensation pressure is assumed 1.4 MPa. A slight subcooling of 0.1 K is set considering the moving boundary of different phase district in transient operation. The mass flow of cooling water is calculated to be 2.26 kg/s by iterative computation to reach the highest system efficiency. Table 1 Main parameters for initial design state.
Parameters
value
Temperature of flue gas
423 K
Mass flow of flue gas
1 kg/s
Journal Pre-proof Temperature of water
298 K
Ambient pressure
101 kPa
Ambient temperature
298 K
Pinch temperature difference in the heater
10 K
Pinch temperature difference in the condenser
8K
High pressure
7000 kPa
Low pressure
1400 kPa
Turbine inlet temperature
413 K
3. Dynamic models The dynamic model of supercritical ORC is developed here to research its dynamic performance. The heater and condenser are represented as typical counter flow straight pipes for the simplicity[17]. In addition, the quasi-stationary models of pump and turbine are employed due to much faster response speed of pump and turbine [15, 17]. The relaxation time of the pump and turbine is so short that they are supposed to operate in a steady state in every moment. The dynamic model of the system is finally created according to the coupling relationships of component modules as shown in Fig.2. Therefore, the equations based on mass and energy conservation can be simultaneous to get a unique solution. The equations are implemented under the Matlab/Simulink platform and the thermo-physical properties of fluids are derived from REFPROP [24].
Journal Pre-proof
Fig.2 Coupling relationships of parameters in each component module
3.1. Heater The Shell-and-tube heat changer is selected for the heater. The working fluid flows inside the tube and the flue gas flows on the shell in the counter-flow heater. The geometric parameters are listed in Table 2. The liquid fluid is heated into vapor in the supercritical condition. It should be noted that the thermal properties of working fluid change sharply with temperature at near-critical state as shown in Fig.3, which are obtained from REFPROP. The change of working fluid density at below 7 MPa is more severe. In addition, at above 8 MPa, the density of working fluid is nearly linear relation with the temperature and there is no sharp decrease in working fluid density. Direct measurements and calculations in the critical region are difficult due to the sharp fluctuations of thermophysical properties of working fluid [25]. Besides, at the
Journal Pre-proof critical point, the classical equations of state fail to describe the asymptotic behavior of thermodynamic properties [26]. The crossover equations of state (EOS) can be used to account for the non-analytical singular behavior in the critical region [27-29]. Many efforts have been made in the development of crossover EOS for mixtures [25, 30, 31]. Oyewunmi [32] employed the statistical associating fluid theory (SAFT) for the non-ideal, asymmetric mixture of ORC system to predict the thermodynamic state properties. However, there is no crossover multiparameter EOS for R134a, R32 and their mixtures so far. The calculation approaches of REFPROP [24] are employed in this work in order to ensure the consistency of the thermal properties involved in the cycle calculation. Considering the reliability and accuracy of the calculation of thermal properties, the thermodynamic states of ORC in this work avoid the tiny regions that are very close to the critical points of R32 and R134a. The high pressure of system in the real-time simulation is higher than 6.7 MPa, the ratio of simulation pressure and the critical pressure of the mixture is higher than 1.3. The fluctuations of thermal properties are much smaller as shown in Fig.3. The FV model cannot realize the identification of near-critical regions and thus cannot ensure the accuracy in the calculations of pseudo-critical region based on dynamic operations. Therefore, an improved dynamic model for zone switching based on density distribution is proposed to achieve high model accuracy in a broad range of operating conditions. The fluctuations of thermal properties in each region are approximately same and the effect of sharp fluctuations of critical region can be
Journal Pre-proof weakened. The schematic diagram of the improved dynamic model is shown in Fig.4. For cell i, i , hi and Ti are the lumped parameters of the cell. Li is the heat transfer length, which varies to ensure uniform thermo-physical properties of working fluids in each cell. The number of discretized cells N is variable, determined by density difference of working fluid in the heater and that in a cell. Fig.5 indicates how the 0 t 1 heater is discretized. In the figure, out is the initial designed value, out is the
density of working fluid at the outlet of heater at the previous moment. The density difference of fluid in each cell ∆ρ is a constant, which is set 50 kg/m3 after weighing the accuracy and computational efficiency. The temperature distribution in cells is not uniform as shown in Fig.4 (b). The temperature difference in the cell around the pseudo-critical point is smaller to ensure accuracy, and thus the variation effects around the pseudo-critical point can be avoided.
Fig.3 (a) Properties of R134a/R32[0.7/0.3] varying with temperature at 5.2MPa and (b) Density of R134a/R32[0.7/0.3] under different temperatures
Journal Pre-proof
Fig.4 (a) Schematic diagram of the improved dynamic model of the heater and (b) the distribution of density and temperature. Table 2 Main geometric parameters of heater and condenser Parameter
Heater
Condenser
Specific heat of tube wall material Cw
470
470
Density of tube wall material ρw (kg/m3)
7790
7790
Tube pitch (m)
0.025
0.025
Diameter of tubes Do (m)
0.019
0.019
Thickness of tubes δ (m)
0.002
0.002
Pipe arrangement
Regular triangle
Regular triangle
(J∙kg-1∙K-1)
Journal Pre-proof
Fig.5 Calculating flow chart of the heater model
The following assumptions for the model are considered to simplify calculation: a) A one-dimensional flow of working fluid in the heat exchange tube is simplified. b) The pressure is assumed to be a constant and thus the momentum balance is ignored. c) The axial heat conduction of fluids and tubes is neglected. d) The effects of gravity potential energy and kinetic energy are not considered. For each discretized zone, the energy and mass conservation equations are established.
Journal Pre-proof Mass conservation of working fluid: Ai m 0 t z
(2)
Energy conservation of working fluid:
Ai h Ai P mh i Di Tw T t z
(3)
Energy conservation of the tube wall: Cw w Aw
Tw i Di Tw T o Do Ths Tw t
(4)
where Af and Aw are the cross section areas of inner pipe and pipe wall, respectively.
i and o is the heat transfer coefficient inside and outside the tube, respectively, which can be calculated [33].
Nu
f / 8 Re 1000 Pr D 1 12.7 f / 8 0.5 Pr 2/3 1
(5)
Integrate the conservation equations (1)-(3) along the axial coordinate for each segment to establish the improved MB dynamic model of heater and the detail results are shown in the Appendix A. The boundary of each cell is moving for dynamic operation conditions, which is determined by the density distribution of working fluid. The lumped value of parameters in each cell is determined by the pressure and the mean of specific enthalpy of fluid at the inlet and outlet of cell, respectively.
3.2. Condenser The Shell-and-tube heat exchanger is selected for the condenser and the main parameters are listed in Table 2. The MB model is employed for the condenser, which is established from three zones characterized by different phases as shown in Fig.6. In
Journal Pre-proof addition, some simplifying assumptions same as that of the heater model are considered.
Fig.6 Schematic diagram of the MB model of the condenser
The moving boundaries are determined by the saturated parameters of real-time pressure. The conservation equations of each region are established as follows the detailed results can be found in Appendix B: Equations of MB model of the condenser. Mass conservation of working fluid: Ai m 0 t z
(6)
Energy conservation of working fluid:
Ai h Ai P mh i Di Twi T t z
(7)
Energy conservation of the tube wall: Cw w Aw
Tw i Di Tw T o Do Tcs Tw t
(8)
In the subcooled and superheated region, the lumped parameters are determined by condensation pressure and the mean of specific enthalpy of fluid at the inlet and
Journal Pre-proof outlet of the region. While for the condensation region, the lumped parameters are calculated as follows [34]:
v 1 l
(9)
h v hv 1 l hl
(10)
where is the void fraction, which is expressed as[35]:
1
(11)
1 1+ 1 l x v
where x is the vapor quality. In addition, the condensation heat transfer coefficient of the zeotropic mixtures side is calculated by equation (8) according to Bell and Ghaly [36] method. 1
1
mix
0.0058 0.557 Pr Y 3.8 l 0.8 0.4 l v 0.023 Rel Prl 1 0.95 D Z 14 v v
Yv xC pv
Tglide
(12)
(13)
hvap
where Tglide is the temperature glide; hvap is the heat of vaporization. Z is Shah's correlating parameter [37]. Z 1/ x 1
0.8
Pr
0.4
(14)
v is the superficial heat transfer coefficient of vapor phase [38]. GxD v 0.023 v where G is the mass velocity.
0.8
0.4 Prv v
D
(15)
Journal Pre-proof 3.3. Pump and turbine The quasi-stationary hypotheses are considered for the pump and turbine. The mass flow of working fluid in pump can be obtained: p vp 3Vcyl p m
(16)
where vp is the volumetric efficiency; Vcyl and p are the cylinder volume and pump speed. Assume that the isentropic efficiency of pump is 0.8 and the specific enthalpy at the outlet of turbine is computed by:
h4
v3 P
p
h3
(17)
where P is the pressure difference of working fluid between the inlet and outlet of the pump. The turbine is simplified as a nozzle [34].
mt ct 1 Ph Pc
(18)
h2 h1 s h1 h2 s
(19)
where cexp is a coefficient; Ph and Pc are the pressure in heater and condensation pressure of working fluid; s is the isentropic efficiency of the turbine assumed to be 0.8. The network of system can be obtained from the equations (17) and (18):
W mt h1 h2 The thermal efficiency of cycle can be calculated by:
(20)
Journal Pre-proof
h1 h2 h1 h4
(21)
The model parameters of pump and turbine are summarized in Table 3. Table 3 Model parameters of pump and turbine
Parameter
Value
Parameter
Value
ηvp
0.6
Vcyl
1.3e-5 m3
ωp
45 rps
ηp
0.8
ct
3.176e-4 m2
𝜂s
0.8
3.4. Model validation A validation of the model in this work has been performed with the experimental ones in Ref [39]. The same operation and heat transfer conditions are assumed. The results when the mass flow rate of heat source and cooling water change same as the literature are compared respectively in Fig.7. The model in this work, as an extension of MB model applied for ORC, can be validated. It reveals that the results of this model agree well with the literature [39]. The small discrepancies occur some of time due to the large uncertainty on the flow rate measurement for the hot water and cooling water, as well as due to the inevitable ambient losses [39].
Journal Pre-proof
Fig.7 Comparison of results between this work and the literature: (a) outlet temperature of heat source as the heat source mass flow rate change; (b) outlet temperature of cooling water as the cooling water mass flow rate change.
4. Results and discussion 4.1. Stability analysis Based on the given boundary conditions, the simulation results of pressures and mass flow and their convergence status are shown in Fig.8. The model of system is converged to a steady state and it has strong robustness. The dynamic characteristics of the system are investigated by introducing disturbances of heat source. The step disturbances of temperature and mass flow of flue gas are tested in this work. In order to compare the response speed of dynamic parameters, a ‘response time’ τ is defined,
Journal Pre-proof which is the minimum time required for the variable to reach the 5% error of the steady-state value and remain within this ranges [40].
Fig.8 Steady state convergence of the system
4.2. Dynamic behavior of the main parameters on variable heat source conditions 4.2.1. Effect of heat source temperature on the parameters of fluid in the heater
A 10K step increase of the heat source temperature is introduced at 5000s when the system has reached a steady state. Fig.9 indicates the dynamic response of the pressure in the heater (Ph). A gradual increase trend is present in the pressure of the heater. The pressure is a key parameter in the heater model, because the lumped parameters of all cells are closely related with it. As a parameter that runs throughout the entire heater from the inlet to the outlet, the response speed of Ph can approximate
Journal Pre-proof the response speed of the heater. The response time of Ph is 1421s, which indicates that the heater of the supercritical ORC has a large thermal inertia. Long time should be taken to settle down when heat source temperature changes.
Fig.9 Dynamic response of the pressure in heater
The mass flow rate of working fluid decreases firstly and then increases as shown in Fig.10 (b). The working flow rate is determined by the turbine as shown in Equation (17), which is related to the fluid density at the inlet of turbine and pressure difference. The pressure ratio and fluid temperature at the inlet of turbine increase along with the increase of Ph (Fig.9). The fluid density of working fluid varies unevenly with temperature due to the effects around the pseudo-critical point as shown in Fig.3. Under the combined influence of Ph and fluid temperature at the turbine inlet, the density first decreases and then increases as shown in Fig.10 (a), and thus a similar trend of the response of mass flow is present.
Journal Pre-proof
Fig.10 Variation of (a) density at the inlet of turbine and (b) mass flow rate of supercritical ORC
Fig.11 illustrates the variation of the temperature of heat source at the heater outlet. The outlet temperature of heat source increases first and then decreases. The heat absorption of working fluid and outlet temperature of heat source vary more slowly than inlet temperature of heat source due to the thermal inertia of the heater. After a period of time, the heat source and the working fluid in the heater have reached a state with a better heat exchange, and the heat absorption of the working medium is increased, so that the heat source outlet temperature decreases. The heat absorption of working fluid increases but the mass flow decreases first and then increases, so the specific enthalpy enhancement of working fluid increases firstly and decreases secondly. Therefore, the specific enthalpy of working fluid at the heater outlet experiences such trend as shown in Fig.12.
Journal Pre-proof
Fig.11 Variation of outlet temperature of heat source of the heater
Fig.12 Specific enthalpy of working fluid at the heater outlet
The thermal efficiency and network increase by 2.44kW and 0.76% as shown in Fig.13. There is a slight decrease in thermal efficiency affected by the fluid specific enthalpy at the heater outlet shown in Fig.12, but it does not occur in the network probably because the variation of mass flow (Fig.10) neutralizes the variation effects of the outlet specific enthalpy of working fluid.
Journal Pre-proof
Fig.13 Variation of (a) thermal efficiency and (b) network of supercritical ORC
Fig.14 indicates the variation of the exergy efficiency of the heater. There is a sudden drop at the beginning. It is because there is a lag in the outlet parameters response of the heater. When a step disturbance is introduced to the inlet temperature of heat source, the outlet parameter does not respond simultaneously and changes in a lag, so the exergy difference of heat source in the heater increases suddenly at first, and thus the exergy efficiency suddenly drops and then increases gradually.
Fig.14 Variation of exergy efficiency of the heater
The system was steady at 4500s, and reached a steady state again at 15000s after introducing a step increase in heat source temperature. Fig.15 shows the proportion of the heat transfer area of each cell in the total area at 4500s and 15000s. The heat
Journal Pre-proof transfer area required can be used as a reference for heat transfer capacity. A lager heat transfer area of a cell should be required if the heat exchange amount in the cell is big or heat transfer coefficient is low. Accordingly, the abnormal situation and the position where is abnormal can be detected. It can be seen that lager heat transfer areas of a cell are mainly distributed in the middle of the heater, which is more obvious at 4500s. Fig.16 indicates the temperature difference of working fluid between the inlet and outlet of a cell along the heater in the flow direction of working fluid. The temperature of working fluid increases in the flow direction, heated by heat source. The thermo-physical properties of the fluids in supercritical pressures are highly variable with temperature as shown in Fig.3, so the temperature difference of a cell around pseudo-critical point is smaller as shown in Fig.16. The location of pseudo-critical point occurs at 77.9% of total heater along the flow direction of working fluid at 4500s and at 69.8% of total heater at 15000s, where should be paid more attention.
Fig.15 Proportion of heat transfer area of each cell in the total area at 4500s and 15000s
Journal Pre-proof
Fig.16 Temperature of working fluid and temperature difference between the inlet and outlet of each cell along the heater (a) t=4500s, (b) t=15000s
4.2.2. Effect of mass flow rate of heat source
Fig.17 describes the variation of Ph, the outlet temperature of heat source and the mass flow rate of working fluid for 10% step increase of exhaust mass flow rate. It can be seen that the effects in the change of mass flow and inlet temperature of heat source are similar. The discrepancies are less than 10%. The effect of heat source temperature can be investigated instead of that of heat source mass flow.
Journal Pre-proof
Fig.17 (a) the pressure in heater Ph, (b) the outlet temperature of heat source and (c) mass flow rate of working fluid in supercritical ORC when the temperature and mass flow of heat source change.
4.3. Influence of the heat transfer coefficients and heat source fluctuations on dynamic response of the system 4.3.1. Relationship between the response time and heat transfer coefficients
It is found that the higher efficiency and the net output power are achieved in supercritical ORC, whereas a lager thermal inertia is exhibited as well. The heater is the main contributor to the thermal inertia of the system, and the heat transfer in the
Journal Pre-proof heater is one of impact factors. The heat transfer coefficients of heat exchangers at the initial design state are listed in Table 4. The heat transfer coefficient in the condenser is much larger than that in the heater due to good heat transfer performance of cold water. The heater in ORC is designed with smooth tubes and the heat transfer coefficient outside the tube is only about 37 W/m2∙K due to the poor heat transfer with gas [41]. Improving heat transfer outside the tubes can effectively enhance the heat exchange of working fluid and the heat source. In order to investigate the effects of heat transfer coefficients in heater on the dynamic response of the system, seven heaters with different heat transfer coefficients outside tubes are tested. 0 is the heat transfer coefficients outside tubes of the original heater in the system I. The heat transfer coefficients outside tubes in the system II-VII is set to 1.5, 2, 2.5, 3, 3.5 and 4 times of that in original heater, respectively, as listed in Table 5. The systems I-VII operate in the same conditions and their network and thermal efficiency are nearly same. The response time of Ph and area of the heater are shown in Fig.18. As the heat transfer coefficient increases, the response time and the required heat transfer area decrease significantly, especially in the case of the initial increase of the heat transfer coefficients. It can be seen that the heat transfer coefficients and response time are approximately exponential. A function of the average of overall heat transfer coefficient in heater and the response time of Ph is fitted as shown in Fig.19. The average of overall heat transfer coefficient is expressed as: U
1 1/ 1/ i 1/ o L L
(22)
Journal Pre-proof where L is the total length of heater, the heat conduction of tube wall and fouling resistance are ignored. The fitting function is obtained by second-order exponential fitting and the confidence coefficient is 95%.
4304 exp(0.04003U ) 439.3exp(0.002232U )
(23)
Table 4 Heat transfer coefficients of heat exchangers at the initial design state
Parameters
Value (W/m2⋅K)
Average αi of the heater
1773.3
Average αo of the heater
37.441
Average αi of the subcooled region in the condenser
1231.2
Average αi of the condensation region in the condenser
7592.2
Average αi of the superheated region in the condenser
909.35
Average αo of the condenser
2551.2
Table 5 Systems with different heat transfer coefficients
System
I
II
III
IV
V
VI
VII
αo
α0
1.5α0
2α0
2.5α0
3α0
3.5α0
4α0
Journal Pre-proof
Fig.18 Heater area and response time of the pressure in heater
Fig.19 Response time of Ph and the prediction function of response time and heat transfer coefficients
The overall heat transfer coefficient of heater can achieve 1286.8 W/m2∙K, when the flue heat source is replaced by hot water. In this case, the response time of Ph in simulation is 31 s, which calculated by the fitting function is 24.9 s. The prediction function turns out to be acceptable. There is a small deviation since the fitting
Journal Pre-proof function is predicted under the exhaust heat source conditions.
4.3.2. Effect of heat source fluctuation on different systems
The system cannot respond instantaneously to changes in heat source parameters because the thermal inertia and the response time depend largely on the heat transfer characteristics. Therefore, some sudden fluctuation peak in high frequency of heat source can be filtered. Depending on the response time to the different frequency input signal, a “dynamic regime number” Γ [40] is defined as:
1/ (2 fsinusoidal )
(24)
where fsinusoidal is the oscillation time of a sinusoidal profile. Three different “dynamic regimes” can be categorized due to different dynamic regime numbers [40] as shown in Fig.20.
Fig.20 Definition of dynamic regimes
The dynamic behavior of Ph in three systems with different thermal inertia is shown in Fig.21, in the case that the fluctuation of heat source temperature is sinusoidal with a frequency of 0.001 Hz and amplitude of 10. The sinusoidal fluctuation is applied representatively, since the transformations can be developed by
Journal Pre-proof Fourier analysis for other fluctuation profiles. In comparison, the system I, III and V with heat transfer coefficients outside tubes of 0 , 2 0 and 3 0 mentioned above are selected. The dynamic regime number of the systems I, III and V is 2.8, 1.3 and 0.85, respectively. The responses of system I and III fall into the dynamic regime II and that of system V is in the dynamic regime I due to its smaller inertia. It can be seen that the system I is affected less by the fluctuation of heat source, but the system responds slowly when the demand power changes.
Fig.21 Response of Ph with different heat transfer coefficients
Fig.22 indicates the response of system VII in the cases of the sinusoidal fluctuation of heat source temperature with frequencies of 0.001 Hz, 0.01 Hz and 0.1Hz, respectively. When the frequency of heat source is 0.001 Hz, the response of Ph has a larger fluctuation amplitude and time lag. In contrast, there is almost no
Journal Pre-proof fluctuation of the response as the frequency of heat source is 0.1 Hz. It is indicated that there is little high-frequency heat source fluctuations effect on large inertia systems. It can be concluded that different dynamic behaviors are demonstrated due to the systems with different inertia and the different fluctuations of heat source. Different strategies can be taken on the basis of different demands and application conditions. A larger inertia system could be selected for the requirement of stability, or the heat exchange could be enhanced for a faster response.
Fig.22 Response of Ph when the frequency of heat source fluctuation is 0.001Hz, 0.01Hz, and 0.1Hz.
5. Conclusion In this work, an improved dynamic model based on the density distribution is developed. The dynamic behaviors of supercritical ORC using mixture working fluid
Journal Pre-proof are presented and the comparisons between the performance of supercritical ORC and subcritical ORC are conducted. The main conclusions are as follows: (1) The improved model can realize the identification of near-critical regions and achieve higher accuracy of model in a broad range of operating conditions. In the case of 10 K step change of heat source temperature, some fluctuation may occur in the outlet parameters of the heater and condenser due to the effects around pseudo-critical point and the thermal inertia. The position of pseudo-critical point is detected, which occurs at 77.9% of total heater along the flow direction of working fluid at 4500s and at 69.8% of total heater at 15000s. (2) The response time and the area required of the heater can be effectively reduced, by enhancing heat transfer performance of the heater. The response time of the pressure and heat transfer coefficients in the heater are approximately exponential. A prediction function has been fitted, which gives guidelines for a design of supercritical ORC application under the real-time conditions and a control system to improve stability performance. (3) The systems with large thermal inertia can remain stable under the high frequency fluctuations of heat source temperature, and the systems with small thermal inertia will be in quasi-stable state under the low frequency fluctuations of heat source temperature. The response characteristic can be classified three regimes roughly. Generally, the amplitude of sinusoidal
Journal Pre-proof response is a fraction of that of step response, and the amplitude is negatively correlated with the frequency and positively correlated with the system thermal inertia. The method of dynamic analysis in this work can be extended for different applications and control guidelines. In the future, more efforts and attention will be paid out to accurately calculate the thermal properties of working fluid at critical point and improve the accuracy of the model.
Journal Pre-proof Nomenclature
A
area (m2)
Cp
specific heat (J/kg∙K)
c
a coefficient
D
diameter (m)
f
friction factor
G
mass velocity (kg/s)
h
specific enthalpy (kJ/kg)
L
length (m)
m
mass flow rate (kg/s)
Nu
Nuselt number
P
pressure (kPa)
Pr
Prandtl number
Re
Reynolds number
T
temperature (K)
t
time (s)
V
volume (m3)
W
work (kW)
x
vapor quality
Journal Pre-proof Greek letters ρ
density (kg/m3)
α
heat transfer coefficient (W/m2∙K)
γ
void fraction
λ
thermal conductivity (W/m∙K)
μ
dynamic viscosity (Pa∙s)
η
efficiency
τ
response time (s)
Г
dynamic regime
ω
revolution speed (rps)
Subscripts and superscripts amb
Ambience
c
condensation
cs
cooling sink
cyl
cylinder volume
h
heater
f
fluid
glide
glide
hs
heat source
i
inside
in
inlet
Journal Pre-proof l
liquid
mix
mixture
o
outside
out
outlet
p
pump
s
isentropic
sub
subcooled
super
superheat
tp
two-phase
t
turbine
v
vapor
w
wall
1-8
state point
Abbreviation ORC
Organic Rankine Cycle
MB
Moving Boundary
FV
Finite Volume
Acknowledgements This work is supported by National Natural Science Foundation of China (No. 51576019), the “artificial intelligence” key project of Chongqing (No. cstc2017rgzn-zdyfX0030).
Journal Pre-proof
Appendix A: Derivation of the heater model equations. Mass conservation of working fluid: Ai m 0 t z
(A.26)
Energy conservation of working fluid:
Ai h Ai P mh i Di Tw T t z
(A.27)
Energy conservation of the tube wall: Cw w Aw
Tw i Di Tw T o Do Ths Tw t
(A.28)
Integrate Eqs.(A.1)-(A.3) in each region respectively. It is assumed that z is the length that increases in the direction of flow as shown in Fig.4 (a). The Leibniz’s rule (Eq.(A.4)) is applied on the conservation equations to simplify them.
z2
z1
f z , t dz dz d z2 dz f z , t dz f z2 , t 2 f z1 , t 1 t dt z1 dt dt
(A.29)
The lumped parameters of each cell are calculated by: hi hi _ in hi _ out 2 i f hi , P Ti f hi , P
(A.30)
(1) Cell 1 The density of working fluid at the outlet of Cell 1:
1_out in Mass conservation of working fluid:
(A.31)
Journal Pre-proof dL A L1 1 1 1_ out 1 min m1_ out dt t
(A.32)
Energy conservation of working fluid:
dh d dP dL1 A 1 1 h1 1 min hin m1_ out h1_ out i1 Di Tw1 T1 L1 1h1 1_ out h1_ out dt dt dt dt (A.33) Energy conservation of the tube wall: Cw w Aw
dTw1 i1 Di T1 Tw1 o1 Do Ths1 Tw1 dt
(A.34)
(2) Cell i (1
i_o ut i _ in in i
(A.35)
Mass conservation of working fluid: d Li i dL A Li i _ in 1_ out i 1 i i_ out i t dt dt
m m i _ in i_ out
(A.36)
Energy conservation of working fluid:
d i dhi i dP A i hi hi 1 Li dt dt P dt d Li dL A i _ in hi _ in i_ out hi_ out i 1 A i hi i _ out hi _ out i dt dt mi _ in hi _ in mi_ out hi_ out ii Di Twi Ti
(A.37)
Energy conservation of the tube wall: Cw w Aw
dTwi ii Di Ti Twi o1 Do Thsi Twi dt
(3) Cell N The heat transfer length of Cell N:
(A.38)
Journal Pre-proof LN L Li
(A.39)
N 1
Mass conservation of working fluid: d Li A N _ in N N 1 LN N dt P mN _ in mout
1 N 2 h h
P
dhN _ in dP 1 LN N dP dt 2 h
P
dhout dt
(A.40) Energy conservation of working fluid:
A N _ in hN _ in N hN
d Li N 1
dt
N 1 ALN N hN 2 hout
dhout dt P dP hN N 1 P h dt
1 dh dhN _ in 1 ALN N N _ in hN N 2 dP 2 h P dP mN _ in hN _ in mout hout iN Di LN TwN TN
(A.41)
Energy conservation of the tube wall: Cw w Aw
dTwN iN Di TN TwN oN Do ThsN TwN dt
(A.42)
Appendix B: Equations of MB model of the condenser. (1). Superheated region: Eqs.(A,1)-(A.3) are integrated over the subcooled region from 0 to L0 as shown in Fig.6 and are simplified with Leibniz’s rule. The lumped parameters of superheated region are calculated by:
hsuper hin hv 2 super f hsuper , P Tsuper f hsuper , P
(A.43)
Journal Pre-proof Mass conservation of working fluid:
dL 1 super A super v 0 L0 super P dt 2 hv h
P
dhv dP 1 super L0 dP dt 2 hin
P
dhin min m1 dt (A.44)
Energy conservation of working fluid: A super hsuper v hv
dL0 1 A L0 super hsuper super dt 2 hin
1 dh 1 AL0 super v hsuper super 2 dP 2 hv
P
dhv hsuper dP
min hin m1hv i 0 Di L0 Tw0 Tsuper
dhin dt P dP super 1 P h dt
(A.45)
Energy conservation of the wall:
dT T T dL Cw w Aw w0 w0 w1 0 i 0 Di T0 Tw0 o 0 Do Ths 0 Tw0 L0 dt dt
( A.46)
(2). Two-phase region: Eqs.(A,1)-(A.3) are integrated from L0 to L0 L1 . The lumped parameters of two-phase region are calculated by:
tp v 1 l tp htp v hv 1 l hl
(A.47)
Mass conversation of working fluid:
dL d dP dL d v A v l 0 v l 1 L1 1 l m1 m2 dt dt dP dt dP (A.48) Energy balance equation of working fluid: d l hl dP d v hv dL dL A l hl v hv 0 1 l hl v hv 1 L1 1 1 dt dt dP dP dt
m1hv m2 hl i1 Di L1 Tw1 T1
Journal Pre-proof (A.49) Energy balance equation of the wall: Cw w Aw
dTw1 i1 Di T1 Tw1 o1 Do Ths1 Tw1 dt
(A.50)
(3) Subcooled region: Eqs.(A,1)-(A.3) are integrated from L0 L1 to L. The lumped parameters of superheated region are calculated by:
hsub hout hl 2 sub f hsub , P Tsub f hsub , P
(A.51)
Mass balance equation of working fluid:
dL dL 1 sub A l sub 0 l sub 1 L2 sub P h 2 hl dt dt m2 mout
P
dhl dP 1 L2 dP dt 2 hl
P
dhl dt
(A.52) Energy balance equation of working fluid:
A l hl sub hsub
d L0 L1 1 AL2 sub hsub sub dt 2 hout
dhout dt P
1 dP dh 1 dhl AL2 sub l hsub sub hsub sub 1 dP 2 h '' P dP P h dt 2 m2 hl mout hout i 2 Di L2 Tw 2 T2
(A.53)
Energy balance equation of the wall:
dT T T dL dL Cw w Aw w 2 w1 w 2 1 2 i 2 Di T2 Tw 2 o 2 Do Ths 2 Tw 2 L2 dt dt dt (A.54)
Reference [1] Li J, Ge Z, Duan Y, Yang Z. Design and performance analyses for a novel
Journal Pre-proof organic Rankine cycle with supercritical-subcritical heat absorption process coupling. Applied Energy. 2019;235:1400-14. [2] Vetter C, Wiemer H-J, Kuhn D. Comparison of sub- and supercritical Organic Rankine Cycles for power generation from low-temperature/low-enthalpy geothermal wells, considering specific net power output and efficiency. Applied Thermal Engineering. 2013;51(1-2):871-9. [3] Yagli H, Koc Y, Koc A, Gorgulu A, Tandiroglu A. Parametric optimization and exergetic analysis comparison of subcritical and supercritical organic Rankine cycle (ORC) for biogas fuelled combined heat and power (CHP) engine exhaust gas waste heat. Energy. 2016;111:923-32. [4] Schuster A, Karellas S, Aumann R. Efficiency optimization potential in supercritical Organic Rankine Cycles. Energy. 2010;35(2):1033-9. [5] Dai X, Shi L, An Q, Qian W. Screening of hydrocarbons as supercritical ORCs working fluids by thermal stability. Energy Conversion and Management. 2016;126:632-7. [6] Chen Y, Lundqvist P, Johansson A, Platell P. A comparative study of the carbon dioxide transcritical power cycle compared with an organic rankine cycle with R123 as working fluid in waste heat recovery. Applied Thermal Engineering. 2006;26(17):2142-7. [7] Mikielewicz D, Mikielewicz J. A thermodynamic criterion for selection of working fluid for subcritical and supercritical domestic micro CHP. Applied Thermal Engineering. 2010;30(16):2357-62. [8] Xu H, Gao N, Zhu T. Investigation on the fluid selection and evaporation parametric optimization for sub- and supercritical organic Rankine cycle. Energy. 2016;96:59-68. [9] Chen H, Goswami DY, Rahman MM, Stefanakos EK. A supercritical Rankine cycle using zeotropic mixture working fluids for the conversion of low-grade heat into power. Energy. 2011;36(1):549-55. [10] Radulovic J, Castaneda NIB. On the potential of zeotropic mixtures in supercritical ORC powered by geothermal energy source. Energy Conversion and Management. 2014;88:365-71. [11] Dong B, Xu G, Cai Y, Li H. Analysis of zeotropic mixtures used in high-temperature Organic Rankine cycle. Energy Conversion and Management. 2014;84:253-60. [12] Xi H, Li M-J, He Y-L, Zhang Y-W. Economical evaluation and optimization of organic Rankine cycle with mixture working fluids using R245fa as flame retardant. Applied Thermal Engineering. 2017;113:1056-70. [13] Xie H, Yang C. Dynamic behavior of Rankine cycle system for waste heat recovery of heavy duty diesel engines under driving cycle. Applied Energy. 2013;112:130-41. [14] Quoilin S, Aumann R, Grill A, Schuster A, Lemort V, Spliethoff H. Dynamic modeling and optimal control strategy of waste heat recovery Organic
Journal Pre-proof Rankine Cycles. Applied Energy. 2011;88(6):2183-90. [15] Horst TA, Tegethoff W, Eilts P, Koehler J. Prediction of dynamic Rankine Cycle waste heat recovery performance and fuel saving potential in passenger car applications considering interactions with vehicles’ energy management. Energy Conversion and Management. 2014;78:438-51. [16] Zhang J, Zhang W, Hou G, Fang F. Dynamic modeling and multivariable control of organic Rankine cycles in waste heat utilizing processes. Computers & Mathematics with Applications. 2012;64(5):908-21. [17] 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 Conversion and Management. 2017;148:724-36. [18] Chowdhury JI, Bao Kha N, Thornhill D. Dynamic model of supercritical organic rankine cycle waste heat recovery system for internal combustion engine. International Journal of Automotive Technology. 2017;18(4):589-601. [19] Chowdhury JI, Bao Kha N, Thornhill D, Hu Y, Soulatiantork P, Balta-Ozkan N, et al. Fuzzy Nonlinear Dynamic Evaporator Model in Supercritical Organic Rankine Cycle Waste Heat Recovery Systems. Energies. 2018;11(4). [20] Wang T, Gao N, Zhu T. Investigation on the optimal condensation temperature of supercritical organic Rankine cycle systems considering meteorological parameters. Energy Conversion and Management. 2018;174:54-64. [21] Olumayegun O, Wang M. Dynamic modelling and control of supercritical CO2 power cycle using waste heat from industrial processes. Fuel. 2019;249:89-102. [22] He C, Liu C, Gao H, Xie H, Li Y, Wu S, et al. The optimal evaporation temperature and working fluids for subcritical organic Rankine cycle. Energy. 2012;38(1):136-43. [23] Jian L, Qiang L, Zhong G, Duan Y, Zhen Y. Thermodynamic performance analyses and optimization of subcritical and transcritical organic Rankine cycles using R1234ze(E) for 100–200 °C heat sources. Energy Conversion & Management. 2017;149:140-54. [24] Lemmon EW, Huber ML, Mclinden MO. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 9.1. NIST NSRDS -. 2010. [25] Xu XH, Duan YY, Yang Z. Prediction of the critical properties of binary alkanol + alkane mixtures using a crossover CPA equation of state. Fluid Phase Equilibria. 2011;309(2):168-73. [26] Yang F, Liu Q, Duan Y, Yang Z. Crossover multiparameter equation of state: General procedure and demonstration with carbon dioxide. Fluid Phase Equilibria. 2019;494:161-71. [27] Kiselev SB. Cubic crossover equation of state1Contribution of the National Institute of Standards and Technology, not subject to copyright in the United States, November 26, 1997.1. Fluid Phase Equilibria. 1998;147(1):7-23. [28] Xu XH, Duan YY. Crossover CPA equation of state for associating fluids.
Journal Pre-proof Fluid Phase Equilibria. 2010;290(1):148-52. [29] Dicko M, Coquelet C. Application of a new crossover treatment to a generalized cubic equation of state. Fluid Phase Equilibria. 2011;302(1):241-8. [30] Shen A-J, Liu Q, Duan Y-Y, Yang Z. Crossover VTSRK equation of state for selected alkane + alkane and CO2 + alkane binary mixtures. Fluid Phase Equilibria. 2016;408:180-9. [31] Belyakov MY, Kulikov VD, Muratov AR, Voronov VP. Crossover Equation of State of a Multi-Component Fluid Mixture in The Vicinity of Liquid-Vapor Critical Points. Chemical Physics. 2018;513:S0301010418303173. [32] Oyewunmi OA, Taleb AI, Haslam AJ, Markides C. An assessment of working-fluid mixtures using SAFT-VR Mie for use in organic rankine cycle systems for waste-heat recovery. 2014;6(4):301-16. [33] Gnielinski V. New equations for heat and mass transfer in turbulent pipe and channel flow. International Chemical Engineering. 1976;16:359-68. [34] Deutsches MO, Elmqvist H, Ab D, Otter M, Jaschinski A, Schweiger C, et al. Moving Boundary Models for Dynamic Simulations of Two-phase Flows. The second international modelica conference. Germany, March 18-19; 2002. [35] Zivi SM. Estimation of Steady-State Steam Void-Fraction by Means of the Principle of Minimum Entropy Production. Transasme Serc. 1964;86(2):247. [36] Bell KJ, Ghaly MA. An Approximate Generalized Design Method for Multicomponent/Partial Condensers. AIChE Symp Ser. 1973;69:72-9. [37] Shah MM. An Improved and Extended General Correlation for Heat Transfer During Condensation in Plain Tubes. HVAC&R Research. 2009;15(5):889-913. [38] Le VL, Kheiri A, Feidt M, Pelloux-Prayer S. Thermodynamic and economic optimizations of a waste heat to power plant driven by a subcritical ORC (Organic Rankine Cycle) using pure or zeotropic working fluid. Energy. 2014;78:622-38. [39] Zhang Y, Deng S, Zhao L, Lin S, Bai M, Wang W, et al. Dynamic test and verification of model-guided ORC system. Energy Conversion and Management. 2019;186:349-67. [40] Jiménez-Arreola M, Pili R, Wieland C, Romagnoli A. Analysis and comparison of dynamic behavior of heat exchangers for direct evaporation in ORC waste heat recovery applications from fluctuating sources. Applied Energy. 2018;216:724-40. [41] Dalkilic AS, Yildiz S, Wongwises S. Experimental investigation of convective heat transfer coefficient during downward laminar flow condensation of R134a in a vertical smooth tube. International Journal of Heat and Mass Transfer. 2009;52(1):142-50.
Journal Pre-proof Highlights
The improved dynamic model with higher accuracy and applicability is developed. The effects around pseudo-critical point of working fluid are investigated. A prediction function of response time and heat transfer coefficients is found. The effects of heat source frequencies and system thermal inertia are researched.
Xiaoxue Chen, Chao Liu, Qibin Li, Xurong Wang, Shukun Wang PII:
S0360-5442(19)32271-6
DOI:
https://doi.org/10.1016/j.energy.2019.116576
Reference:
EGY 116576
To appear in:
Energy
Received Date:
06 August 2019
Accepted Date:
17 November 2019
Please cite this article as: Xiaoxue Chen, Chao Liu, Qibin Li, Xurong Wang, Shukun Wang, Dynamic behavior of supercritical organic Rankine cycle using zeotropic mixture working fluids, Energy (2019), https://doi.org/10.1016/j.energy.2019.116576
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.
Journal Pre-proof
Dynamic behavior of supercritical organic Rankine cycle using zeotropic mixture working fluids
Xiaoxue Chen, Chao Liu*, Qibin Li*, Xurong Wang, Shukun Wang
Key laboratory of low-grade Energy Utilization Technologies and Systems, Ministry of Education, School of Energy and Power Engineering, Chongqing University, Chongqing 400030, China *Corresponding author: 1. Chao Liu. E-mail address: [email protected] 2. Qibin Li. E-mail address: [email protected]
Abstract As a prominent technology for recovering low-grade waste heat, supercritical organic Rankine cycle (ORC) exhibits a better performance due to the higher endothermic temperature and better thermal match with the heat source. The irreversibility in system can be reduced greatly when the zeotropic mixtures are used as working fluid to match the heat source and sink profiles. Affected by the fluctuations in waste heat sources, there is a challenge for ORC to recovery waste heat. An improved dynamic model of supercritical ORC using zeotropic mixture R134a/R32 as working fluid is developed and dynamic behaviors of supercritical ORC are analyzed. It is found that an abnormal fluctuation may occur in some parameters due to the effects around pseudo-critical point and thermal inertia of the
Journal Pre-proof heater. A fitting correlation to predict the response time based on different heat transfer coefficients of heater is found. Besides, three dynamic regimes are defined to investigate the effects of heat source frequencies and system thermal inertia on the dynamic response. As the heat source frequency increases or the heat exchange in the heater is enhanced, the fluctuation amplitude of the pressure decreases in the heater. Key words: Supercritical organic Rankine cycle, zeotropic mixture, dynamic characteristic, thermal inertia, dynamic regime
1. Introduction In recent years, the issue of energy shortage and carbon dioxide emission has been increasingly serious. The utilization of renewable energy source has been attracted more attention to hedge against rising costs of fossil fuels and reduction in greenhouse gas emissions. Furthermore, there is a high potential for waste heat recovery in many industries such as cement industries, chemical industries, iron and steel industries. The organic Rankine cycle (ORC) is a superior technology for utilization of low-grade energy such as solar energy, geothermal energy and waste heat, etc. Generally, two types of ORC for different operation conditions are mainly included: subcritical and supercritical ORC. Li et al. [1] also proposed a novel ORC that couples supercritical and subcritical heat absorption processes, which can better enhance the system efficiency. The comparisons between the performance of
Journal Pre-proof supercritical and subcritical ORC have been conducted in some researches [2, 3]. Supercritical approach, by contrast, can achieve a better performance of system due to the higher endothermic temperature and better thermal match with the heat source [4]. Thermal and exergy efficiencies between the subcritical and supercritical ORC for biogas fueled combined heat and power engine system were compared by Yagli et al. [3]. Their results showed that the supercritical ORC can make much network output and a higher thermal and exergy efficiency of system. The selection of working fluid is crucial to the ORC system performance. Carbon dioxide (CO2) with relatively low critical temperatures and pressures is widely used in supercritical ORC. Some organic fluids such as hydrocarbons [5] and refrigerants R123 [6], R134a [7], R245fa, R143a [8] has also been studied as working fluid of supercritical ORC. Since the phase change of pure working fluids has a large temperature difference to result in a bad thermal match with the heat source or heat sink, the zeotropic mixtures emerge great advantages due to its temperature glide. Chen et al. [9] introduced the idea of using zeotropic mixtures as working fluid for supercritical ORC. The comparisons between the conventional ORC and supercritical ORC indicate that the higher thermal efficiency and exergy efficiency of supercritical ORC can be achieved by using a zeotropic mixture R134a/R32. Radulovic et al. [10] made a comparison on the operation performance of supercritical ORC using six zeotropic mixtures and found the thermal efficiency of system using pure working fluid is overall lower than using zeotropic mixtures. Dong et al.[11] studied zeotropic
Journal Pre-proof mixtures of siloxanes as working fluids in high-temperature ORC. The zeotropic mixture working fluids have been indicated to has a great potential in economic profitability [12]. In order to achieve a better performance, many studies have been carried out, especially for the steady state. In fact, affected by some uncertain factors, the parameters of heat source often fluctuate over the time in actual operating conditions. In addition, the operating mode switching may result in an inefficiency, particularly for the exhaust heat recovery of internal combustion engines [13]. It is significant to predict the dynamic behavior of system under the unsteady state. Quoilin et al. [14] proposed a dynamic model of the ORC focusing specifically on the time-varying performance of heat exchangers. Three different control strategies also were tested and compared to optimize the working conditions. Horst et al. [15] presented an approach for predicting dynamic Rankine Cycle performance and fuel saving potential in passenger car. A multivariable control strategy was proposed in order to achieve both transient performance and steady-state energy saving [16]. Wang et al. [17] established a dynamic model to study the dynamic behaviors and part-load performance of the double loop ORC. The effects of the mass flow rate of high temperature cooling water on the system performance are researched as well. Chowdhury et al. [18] presented a Finite Volume (FV) dynamic model of supercritical ORC and analyzed the dynamic behaviors and transient responses of the system based on the transient heat source. Due to the limitation of time-intensive calculation of FV
Journal Pre-proof model, they also developed a fuzzy-based dynamic model for real-time applications [19]. In view of seasonal and meteorological factors, Wang et al. [20] established a series of supercritical ORC thermodynamic models under the design and off-design conditions. Olumayegun et al. [21] proposed a dynamic model of supercritical CO2 power cycle for cement industry waste heat recovery and found it is preferable to allow the temperature of the turbine inlet to vary with the transient heat source. However, most present dynamic studies are focused on the subcritical system and the literatures on the dynamic analysis of supercritical ORC are seldom. In addition, the thermo-physical properties of the fluids in supercritical pressures are highly variable with the temperature and existing dynamic models cannot eliminate the effects of such variation. This leads to a loss of accuracy or computational efficiency of models. Therefore, an improved model to consider the effect of variable thermo-physical property in supercritical fluid is developed in this work. Both the model accuracy and computational efficiency can be satisfied simultaneously. The dynamic behavior of supercritical ORC is analyzed in this work. The effects around pseudo-critical point and thermal inertia of heater are investigated. The model is discretized by moving boundary, thus variable heat transfer area of a discretized cell can be used as a detection of heat transfer capacity, and furthermore the position of pseudo-critical point in the heater can be found. In addition, the dynamic responses of the supercritical system with different thermal inertia and variable frequencies of heat source are presented. A mathematical correlation between the response time and
Journal Pre-proof transfer coefficients of the heater is found, which can give guidance for control and design of the system on the basis of different demands and application conditions.
2. System description and cycle process Xu [8] pointed out that the critical temperature of most suitable working fluid is 40-65 ℃ lower than the flue gas inlet temperature and the optimum condition is supercritical. Taking the fluids’ properties and operation conditions into consideration, the mixture of R134a/R32 (mass fraction is [0.7/0.3]) is selected in this study, which has been turned out to achieve a high thermal and exergy efficiency of system [9]. The schematic layout and T-s diagram of the supercritical ORC with zeotropic mixture of 0.7R134a/0.3R32 are shown in Fig.1. The working fluid is heated by the flue gas and becomes supercritical vapor directly from the liquid in the heater (4-1). The supercritical vapor enters into the turbine and completes the expansion to generate power. Then the superheated vapor from the turbine flows into the condenser where it is cooled down to subcooled liquid by cooling water but non-isothermally (2-3). After that, the liquid is pressurized (3-4) to be supercritical fluid into the heater, and a new cycle begins.
Journal Pre-proof
Fig.1 (a) Schematic of ORC system and (b) T-S diagram of the process
The main parameters for the initial design point are listed in Table 1 [9, 22]. The critical pressure of 0.7R134a/0.3R32 is 5.13 MPa [9]. The high pressure of initial designed cycle is set 7 MPa and the inlet temperature of working fluid in turbine is set 413 K, which is 10 K lower than the temperature of flue gas [9]. The droplets formation can be avoided in the expansion process [23]. The condensation pressure is assumed 1.4 MPa. A slight subcooling of 0.1 K is set considering the moving boundary of different phase district in transient operation. The mass flow of cooling water is calculated to be 2.26 kg/s by iterative computation to reach the highest system efficiency. Table 1 Main parameters for initial design state.
Parameters
value
Temperature of flue gas
423 K
Mass flow of flue gas
1 kg/s
Journal Pre-proof Temperature of water
298 K
Ambient pressure
101 kPa
Ambient temperature
298 K
Pinch temperature difference in the heater
10 K
Pinch temperature difference in the condenser
8K
High pressure
7000 kPa
Low pressure
1400 kPa
Turbine inlet temperature
413 K
3. Dynamic models The dynamic model of supercritical ORC is developed here to research its dynamic performance. The heater and condenser are represented as typical counter flow straight pipes for the simplicity[17]. In addition, the quasi-stationary models of pump and turbine are employed due to much faster response speed of pump and turbine [15, 17]. The relaxation time of the pump and turbine is so short that they are supposed to operate in a steady state in every moment. The dynamic model of the system is finally created according to the coupling relationships of component modules as shown in Fig.2. Therefore, the equations based on mass and energy conservation can be simultaneous to get a unique solution. The equations are implemented under the Matlab/Simulink platform and the thermo-physical properties of fluids are derived from REFPROP [24].
Journal Pre-proof
Fig.2 Coupling relationships of parameters in each component module
3.1. Heater The Shell-and-tube heat changer is selected for the heater. The working fluid flows inside the tube and the flue gas flows on the shell in the counter-flow heater. The geometric parameters are listed in Table 2. The liquid fluid is heated into vapor in the supercritical condition. It should be noted that the thermal properties of working fluid change sharply with temperature at near-critical state as shown in Fig.3, which are obtained from REFPROP. The change of working fluid density at below 7 MPa is more severe. In addition, at above 8 MPa, the density of working fluid is nearly linear relation with the temperature and there is no sharp decrease in working fluid density. Direct measurements and calculations in the critical region are difficult due to the sharp fluctuations of thermophysical properties of working fluid [25]. Besides, at the
Journal Pre-proof critical point, the classical equations of state fail to describe the asymptotic behavior of thermodynamic properties [26]. The crossover equations of state (EOS) can be used to account for the non-analytical singular behavior in the critical region [27-29]. Many efforts have been made in the development of crossover EOS for mixtures [25, 30, 31]. Oyewunmi [32] employed the statistical associating fluid theory (SAFT) for the non-ideal, asymmetric mixture of ORC system to predict the thermodynamic state properties. However, there is no crossover multiparameter EOS for R134a, R32 and their mixtures so far. The calculation approaches of REFPROP [24] are employed in this work in order to ensure the consistency of the thermal properties involved in the cycle calculation. Considering the reliability and accuracy of the calculation of thermal properties, the thermodynamic states of ORC in this work avoid the tiny regions that are very close to the critical points of R32 and R134a. The high pressure of system in the real-time simulation is higher than 6.7 MPa, the ratio of simulation pressure and the critical pressure of the mixture is higher than 1.3. The fluctuations of thermal properties are much smaller as shown in Fig.3. The FV model cannot realize the identification of near-critical regions and thus cannot ensure the accuracy in the calculations of pseudo-critical region based on dynamic operations. Therefore, an improved dynamic model for zone switching based on density distribution is proposed to achieve high model accuracy in a broad range of operating conditions. The fluctuations of thermal properties in each region are approximately same and the effect of sharp fluctuations of critical region can be
Journal Pre-proof weakened. The schematic diagram of the improved dynamic model is shown in Fig.4. For cell i, i , hi and Ti are the lumped parameters of the cell. Li is the heat transfer length, which varies to ensure uniform thermo-physical properties of working fluids in each cell. The number of discretized cells N is variable, determined by density difference of working fluid in the heater and that in a cell. Fig.5 indicates how the 0 t 1 heater is discretized. In the figure, out is the initial designed value, out is the
density of working fluid at the outlet of heater at the previous moment. The density difference of fluid in each cell ∆ρ is a constant, which is set 50 kg/m3 after weighing the accuracy and computational efficiency. The temperature distribution in cells is not uniform as shown in Fig.4 (b). The temperature difference in the cell around the pseudo-critical point is smaller to ensure accuracy, and thus the variation effects around the pseudo-critical point can be avoided.
Fig.3 (a) Properties of R134a/R32[0.7/0.3] varying with temperature at 5.2MPa and (b) Density of R134a/R32[0.7/0.3] under different temperatures
Journal Pre-proof
Fig.4 (a) Schematic diagram of the improved dynamic model of the heater and (b) the distribution of density and temperature. Table 2 Main geometric parameters of heater and condenser Parameter
Heater
Condenser
Specific heat of tube wall material Cw
470
470
Density of tube wall material ρw (kg/m3)
7790
7790
Tube pitch (m)
0.025
0.025
Diameter of tubes Do (m)
0.019
0.019
Thickness of tubes δ (m)
0.002
0.002
Pipe arrangement
Regular triangle
Regular triangle
(J∙kg-1∙K-1)
Journal Pre-proof
Fig.5 Calculating flow chart of the heater model
The following assumptions for the model are considered to simplify calculation: a) A one-dimensional flow of working fluid in the heat exchange tube is simplified. b) The pressure is assumed to be a constant and thus the momentum balance is ignored. c) The axial heat conduction of fluids and tubes is neglected. d) The effects of gravity potential energy and kinetic energy are not considered. For each discretized zone, the energy and mass conservation equations are established.
Journal Pre-proof Mass conservation of working fluid: Ai m 0 t z
(2)
Energy conservation of working fluid:
Ai h Ai P mh i Di Tw T t z
(3)
Energy conservation of the tube wall: Cw w Aw
Tw i Di Tw T o Do Ths Tw t
(4)
where Af and Aw are the cross section areas of inner pipe and pipe wall, respectively.
i and o is the heat transfer coefficient inside and outside the tube, respectively, which can be calculated [33].
Nu
f / 8 Re 1000 Pr D 1 12.7 f / 8 0.5 Pr 2/3 1
(5)
Integrate the conservation equations (1)-(3) along the axial coordinate for each segment to establish the improved MB dynamic model of heater and the detail results are shown in the Appendix A. The boundary of each cell is moving for dynamic operation conditions, which is determined by the density distribution of working fluid. The lumped value of parameters in each cell is determined by the pressure and the mean of specific enthalpy of fluid at the inlet and outlet of cell, respectively.
3.2. Condenser The Shell-and-tube heat exchanger is selected for the condenser and the main parameters are listed in Table 2. The MB model is employed for the condenser, which is established from three zones characterized by different phases as shown in Fig.6. In
Journal Pre-proof addition, some simplifying assumptions same as that of the heater model are considered.
Fig.6 Schematic diagram of the MB model of the condenser
The moving boundaries are determined by the saturated parameters of real-time pressure. The conservation equations of each region are established as follows the detailed results can be found in Appendix B: Equations of MB model of the condenser. Mass conservation of working fluid: Ai m 0 t z
(6)
Energy conservation of working fluid:
Ai h Ai P mh i Di Twi T t z
(7)
Energy conservation of the tube wall: Cw w Aw
Tw i Di Tw T o Do Tcs Tw t
(8)
In the subcooled and superheated region, the lumped parameters are determined by condensation pressure and the mean of specific enthalpy of fluid at the inlet and
Journal Pre-proof outlet of the region. While for the condensation region, the lumped parameters are calculated as follows [34]:
v 1 l
(9)
h v hv 1 l hl
(10)
where is the void fraction, which is expressed as[35]:
1
(11)
1 1+ 1 l x v
where x is the vapor quality. In addition, the condensation heat transfer coefficient of the zeotropic mixtures side is calculated by equation (8) according to Bell and Ghaly [36] method. 1
1
mix
0.0058 0.557 Pr Y 3.8 l 0.8 0.4 l v 0.023 Rel Prl 1 0.95 D Z 14 v v
Yv xC pv
Tglide
(12)
(13)
hvap
where Tglide is the temperature glide; hvap is the heat of vaporization. Z is Shah's correlating parameter [37]. Z 1/ x 1
0.8
Pr
0.4
(14)
v is the superficial heat transfer coefficient of vapor phase [38]. GxD v 0.023 v where G is the mass velocity.
0.8
0.4 Prv v
D
(15)
Journal Pre-proof 3.3. Pump and turbine The quasi-stationary hypotheses are considered for the pump and turbine. The mass flow of working fluid in pump can be obtained: p vp 3Vcyl p m
(16)
where vp is the volumetric efficiency; Vcyl and p are the cylinder volume and pump speed. Assume that the isentropic efficiency of pump is 0.8 and the specific enthalpy at the outlet of turbine is computed by:
h4
v3 P
p
h3
(17)
where P is the pressure difference of working fluid between the inlet and outlet of the pump. The turbine is simplified as a nozzle [34].
mt ct 1 Ph Pc
(18)
h2 h1 s h1 h2 s
(19)
where cexp is a coefficient; Ph and Pc are the pressure in heater and condensation pressure of working fluid; s is the isentropic efficiency of the turbine assumed to be 0.8. The network of system can be obtained from the equations (17) and (18):
W mt h1 h2 The thermal efficiency of cycle can be calculated by:
(20)
Journal Pre-proof
h1 h2 h1 h4
(21)
The model parameters of pump and turbine are summarized in Table 3. Table 3 Model parameters of pump and turbine
Parameter
Value
Parameter
Value
ηvp
0.6
Vcyl
1.3e-5 m3
ωp
45 rps
ηp
0.8
ct
3.176e-4 m2
𝜂s
0.8
3.4. Model validation A validation of the model in this work has been performed with the experimental ones in Ref [39]. The same operation and heat transfer conditions are assumed. The results when the mass flow rate of heat source and cooling water change same as the literature are compared respectively in Fig.7. The model in this work, as an extension of MB model applied for ORC, can be validated. It reveals that the results of this model agree well with the literature [39]. The small discrepancies occur some of time due to the large uncertainty on the flow rate measurement for the hot water and cooling water, as well as due to the inevitable ambient losses [39].
Journal Pre-proof
Fig.7 Comparison of results between this work and the literature: (a) outlet temperature of heat source as the heat source mass flow rate change; (b) outlet temperature of cooling water as the cooling water mass flow rate change.
4. Results and discussion 4.1. Stability analysis Based on the given boundary conditions, the simulation results of pressures and mass flow and their convergence status are shown in Fig.8. The model of system is converged to a steady state and it has strong robustness. The dynamic characteristics of the system are investigated by introducing disturbances of heat source. The step disturbances of temperature and mass flow of flue gas are tested in this work. In order to compare the response speed of dynamic parameters, a ‘response time’ τ is defined,
Journal Pre-proof which is the minimum time required for the variable to reach the 5% error of the steady-state value and remain within this ranges [40].
Fig.8 Steady state convergence of the system
4.2. Dynamic behavior of the main parameters on variable heat source conditions 4.2.1. Effect of heat source temperature on the parameters of fluid in the heater
A 10K step increase of the heat source temperature is introduced at 5000s when the system has reached a steady state. Fig.9 indicates the dynamic response of the pressure in the heater (Ph). A gradual increase trend is present in the pressure of the heater. The pressure is a key parameter in the heater model, because the lumped parameters of all cells are closely related with it. As a parameter that runs throughout the entire heater from the inlet to the outlet, the response speed of Ph can approximate
Journal Pre-proof the response speed of the heater. The response time of Ph is 1421s, which indicates that the heater of the supercritical ORC has a large thermal inertia. Long time should be taken to settle down when heat source temperature changes.
Fig.9 Dynamic response of the pressure in heater
The mass flow rate of working fluid decreases firstly and then increases as shown in Fig.10 (b). The working flow rate is determined by the turbine as shown in Equation (17), which is related to the fluid density at the inlet of turbine and pressure difference. The pressure ratio and fluid temperature at the inlet of turbine increase along with the increase of Ph (Fig.9). The fluid density of working fluid varies unevenly with temperature due to the effects around the pseudo-critical point as shown in Fig.3. Under the combined influence of Ph and fluid temperature at the turbine inlet, the density first decreases and then increases as shown in Fig.10 (a), and thus a similar trend of the response of mass flow is present.
Journal Pre-proof
Fig.10 Variation of (a) density at the inlet of turbine and (b) mass flow rate of supercritical ORC
Fig.11 illustrates the variation of the temperature of heat source at the heater outlet. The outlet temperature of heat source increases first and then decreases. The heat absorption of working fluid and outlet temperature of heat source vary more slowly than inlet temperature of heat source due to the thermal inertia of the heater. After a period of time, the heat source and the working fluid in the heater have reached a state with a better heat exchange, and the heat absorption of the working medium is increased, so that the heat source outlet temperature decreases. The heat absorption of working fluid increases but the mass flow decreases first and then increases, so the specific enthalpy enhancement of working fluid increases firstly and decreases secondly. Therefore, the specific enthalpy of working fluid at the heater outlet experiences such trend as shown in Fig.12.
Journal Pre-proof
Fig.11 Variation of outlet temperature of heat source of the heater
Fig.12 Specific enthalpy of working fluid at the heater outlet
The thermal efficiency and network increase by 2.44kW and 0.76% as shown in Fig.13. There is a slight decrease in thermal efficiency affected by the fluid specific enthalpy at the heater outlet shown in Fig.12, but it does not occur in the network probably because the variation of mass flow (Fig.10) neutralizes the variation effects of the outlet specific enthalpy of working fluid.
Journal Pre-proof
Fig.13 Variation of (a) thermal efficiency and (b) network of supercritical ORC
Fig.14 indicates the variation of the exergy efficiency of the heater. There is a sudden drop at the beginning. It is because there is a lag in the outlet parameters response of the heater. When a step disturbance is introduced to the inlet temperature of heat source, the outlet parameter does not respond simultaneously and changes in a lag, so the exergy difference of heat source in the heater increases suddenly at first, and thus the exergy efficiency suddenly drops and then increases gradually.
Fig.14 Variation of exergy efficiency of the heater
The system was steady at 4500s, and reached a steady state again at 15000s after introducing a step increase in heat source temperature. Fig.15 shows the proportion of the heat transfer area of each cell in the total area at 4500s and 15000s. The heat
Journal Pre-proof transfer area required can be used as a reference for heat transfer capacity. A lager heat transfer area of a cell should be required if the heat exchange amount in the cell is big or heat transfer coefficient is low. Accordingly, the abnormal situation and the position where is abnormal can be detected. It can be seen that lager heat transfer areas of a cell are mainly distributed in the middle of the heater, which is more obvious at 4500s. Fig.16 indicates the temperature difference of working fluid between the inlet and outlet of a cell along the heater in the flow direction of working fluid. The temperature of working fluid increases in the flow direction, heated by heat source. The thermo-physical properties of the fluids in supercritical pressures are highly variable with temperature as shown in Fig.3, so the temperature difference of a cell around pseudo-critical point is smaller as shown in Fig.16. The location of pseudo-critical point occurs at 77.9% of total heater along the flow direction of working fluid at 4500s and at 69.8% of total heater at 15000s, where should be paid more attention.
Fig.15 Proportion of heat transfer area of each cell in the total area at 4500s and 15000s
Journal Pre-proof
Fig.16 Temperature of working fluid and temperature difference between the inlet and outlet of each cell along the heater (a) t=4500s, (b) t=15000s
4.2.2. Effect of mass flow rate of heat source
Fig.17 describes the variation of Ph, the outlet temperature of heat source and the mass flow rate of working fluid for 10% step increase of exhaust mass flow rate. It can be seen that the effects in the change of mass flow and inlet temperature of heat source are similar. The discrepancies are less than 10%. The effect of heat source temperature can be investigated instead of that of heat source mass flow.
Journal Pre-proof
Fig.17 (a) the pressure in heater Ph, (b) the outlet temperature of heat source and (c) mass flow rate of working fluid in supercritical ORC when the temperature and mass flow of heat source change.
4.3. Influence of the heat transfer coefficients and heat source fluctuations on dynamic response of the system 4.3.1. Relationship between the response time and heat transfer coefficients
It is found that the higher efficiency and the net output power are achieved in supercritical ORC, whereas a lager thermal inertia is exhibited as well. The heater is the main contributor to the thermal inertia of the system, and the heat transfer in the
Journal Pre-proof heater is one of impact factors. The heat transfer coefficients of heat exchangers at the initial design state are listed in Table 4. The heat transfer coefficient in the condenser is much larger than that in the heater due to good heat transfer performance of cold water. The heater in ORC is designed with smooth tubes and the heat transfer coefficient outside the tube is only about 37 W/m2∙K due to the poor heat transfer with gas [41]. Improving heat transfer outside the tubes can effectively enhance the heat exchange of working fluid and the heat source. In order to investigate the effects of heat transfer coefficients in heater on the dynamic response of the system, seven heaters with different heat transfer coefficients outside tubes are tested. 0 is the heat transfer coefficients outside tubes of the original heater in the system I. The heat transfer coefficients outside tubes in the system II-VII is set to 1.5, 2, 2.5, 3, 3.5 and 4 times of that in original heater, respectively, as listed in Table 5. The systems I-VII operate in the same conditions and their network and thermal efficiency are nearly same. The response time of Ph and area of the heater are shown in Fig.18. As the heat transfer coefficient increases, the response time and the required heat transfer area decrease significantly, especially in the case of the initial increase of the heat transfer coefficients. It can be seen that the heat transfer coefficients and response time are approximately exponential. A function of the average of overall heat transfer coefficient in heater and the response time of Ph is fitted as shown in Fig.19. The average of overall heat transfer coefficient is expressed as: U
1 1/ 1/ i 1/ o L L
(22)
Journal Pre-proof where L is the total length of heater, the heat conduction of tube wall and fouling resistance are ignored. The fitting function is obtained by second-order exponential fitting and the confidence coefficient is 95%.
4304 exp(0.04003U ) 439.3exp(0.002232U )
(23)
Table 4 Heat transfer coefficients of heat exchangers at the initial design state
Parameters
Value (W/m2⋅K)
Average αi of the heater
1773.3
Average αo of the heater
37.441
Average αi of the subcooled region in the condenser
1231.2
Average αi of the condensation region in the condenser
7592.2
Average αi of the superheated region in the condenser
909.35
Average αo of the condenser
2551.2
Table 5 Systems with different heat transfer coefficients
System
I
II
III
IV
V
VI
VII
αo
α0
1.5α0
2α0
2.5α0
3α0
3.5α0
4α0
Journal Pre-proof
Fig.18 Heater area and response time of the pressure in heater
Fig.19 Response time of Ph and the prediction function of response time and heat transfer coefficients
The overall heat transfer coefficient of heater can achieve 1286.8 W/m2∙K, when the flue heat source is replaced by hot water. In this case, the response time of Ph in simulation is 31 s, which calculated by the fitting function is 24.9 s. The prediction function turns out to be acceptable. There is a small deviation since the fitting
Journal Pre-proof function is predicted under the exhaust heat source conditions.
4.3.2. Effect of heat source fluctuation on different systems
The system cannot respond instantaneously to changes in heat source parameters because the thermal inertia and the response time depend largely on the heat transfer characteristics. Therefore, some sudden fluctuation peak in high frequency of heat source can be filtered. Depending on the response time to the different frequency input signal, a “dynamic regime number” Γ [40] is defined as:
1/ (2 fsinusoidal )
(24)
where fsinusoidal is the oscillation time of a sinusoidal profile. Three different “dynamic regimes” can be categorized due to different dynamic regime numbers [40] as shown in Fig.20.
Fig.20 Definition of dynamic regimes
The dynamic behavior of Ph in three systems with different thermal inertia is shown in Fig.21, in the case that the fluctuation of heat source temperature is sinusoidal with a frequency of 0.001 Hz and amplitude of 10. The sinusoidal fluctuation is applied representatively, since the transformations can be developed by
Journal Pre-proof Fourier analysis for other fluctuation profiles. In comparison, the system I, III and V with heat transfer coefficients outside tubes of 0 , 2 0 and 3 0 mentioned above are selected. The dynamic regime number of the systems I, III and V is 2.8, 1.3 and 0.85, respectively. The responses of system I and III fall into the dynamic regime II and that of system V is in the dynamic regime I due to its smaller inertia. It can be seen that the system I is affected less by the fluctuation of heat source, but the system responds slowly when the demand power changes.
Fig.21 Response of Ph with different heat transfer coefficients
Fig.22 indicates the response of system VII in the cases of the sinusoidal fluctuation of heat source temperature with frequencies of 0.001 Hz, 0.01 Hz and 0.1Hz, respectively. When the frequency of heat source is 0.001 Hz, the response of Ph has a larger fluctuation amplitude and time lag. In contrast, there is almost no
Journal Pre-proof fluctuation of the response as the frequency of heat source is 0.1 Hz. It is indicated that there is little high-frequency heat source fluctuations effect on large inertia systems. It can be concluded that different dynamic behaviors are demonstrated due to the systems with different inertia and the different fluctuations of heat source. Different strategies can be taken on the basis of different demands and application conditions. A larger inertia system could be selected for the requirement of stability, or the heat exchange could be enhanced for a faster response.
Fig.22 Response of Ph when the frequency of heat source fluctuation is 0.001Hz, 0.01Hz, and 0.1Hz.
5. Conclusion In this work, an improved dynamic model based on the density distribution is developed. The dynamic behaviors of supercritical ORC using mixture working fluid
Journal Pre-proof are presented and the comparisons between the performance of supercritical ORC and subcritical ORC are conducted. The main conclusions are as follows: (1) The improved model can realize the identification of near-critical regions and achieve higher accuracy of model in a broad range of operating conditions. In the case of 10 K step change of heat source temperature, some fluctuation may occur in the outlet parameters of the heater and condenser due to the effects around pseudo-critical point and the thermal inertia. The position of pseudo-critical point is detected, which occurs at 77.9% of total heater along the flow direction of working fluid at 4500s and at 69.8% of total heater at 15000s. (2) The response time and the area required of the heater can be effectively reduced, by enhancing heat transfer performance of the heater. The response time of the pressure and heat transfer coefficients in the heater are approximately exponential. A prediction function has been fitted, which gives guidelines for a design of supercritical ORC application under the real-time conditions and a control system to improve stability performance. (3) The systems with large thermal inertia can remain stable under the high frequency fluctuations of heat source temperature, and the systems with small thermal inertia will be in quasi-stable state under the low frequency fluctuations of heat source temperature. The response characteristic can be classified three regimes roughly. Generally, the amplitude of sinusoidal
Journal Pre-proof response is a fraction of that of step response, and the amplitude is negatively correlated with the frequency and positively correlated with the system thermal inertia. The method of dynamic analysis in this work can be extended for different applications and control guidelines. In the future, more efforts and attention will be paid out to accurately calculate the thermal properties of working fluid at critical point and improve the accuracy of the model.
Journal Pre-proof Nomenclature
A
area (m2)
Cp
specific heat (J/kg∙K)
c
a coefficient
D
diameter (m)
f
friction factor
G
mass velocity (kg/s)
h
specific enthalpy (kJ/kg)
L
length (m)
m
mass flow rate (kg/s)
Nu
Nuselt number
P
pressure (kPa)
Pr
Prandtl number
Re
Reynolds number
T
temperature (K)
t
time (s)
V
volume (m3)
W
work (kW)
x
vapor quality
Journal Pre-proof Greek letters ρ
density (kg/m3)
α
heat transfer coefficient (W/m2∙K)
γ
void fraction
λ
thermal conductivity (W/m∙K)
μ
dynamic viscosity (Pa∙s)
η
efficiency
τ
response time (s)
Г
dynamic regime
ω
revolution speed (rps)
Subscripts and superscripts amb
Ambience
c
condensation
cs
cooling sink
cyl
cylinder volume
h
heater
f
fluid
glide
glide
hs
heat source
i
inside
in
inlet
Journal Pre-proof l
liquid
mix
mixture
o
outside
out
outlet
p
pump
s
isentropic
sub
subcooled
super
superheat
tp
two-phase
t
turbine
v
vapor
w
wall
1-8
state point
Abbreviation ORC
Organic Rankine Cycle
MB
Moving Boundary
FV
Finite Volume
Acknowledgements This work is supported by National Natural Science Foundation of China (No. 51576019), the “artificial intelligence” key project of Chongqing (No. cstc2017rgzn-zdyfX0030).
Journal Pre-proof
Appendix A: Derivation of the heater model equations. Mass conservation of working fluid: Ai m 0 t z
(A.26)
Energy conservation of working fluid:
Ai h Ai P mh i Di Tw T t z
(A.27)
Energy conservation of the tube wall: Cw w Aw
Tw i Di Tw T o Do Ths Tw t
(A.28)
Integrate Eqs.(A.1)-(A.3) in each region respectively. It is assumed that z is the length that increases in the direction of flow as shown in Fig.4 (a). The Leibniz’s rule (Eq.(A.4)) is applied on the conservation equations to simplify them.
z2
z1
f z , t dz dz d z2 dz f z , t dz f z2 , t 2 f z1 , t 1 t dt z1 dt dt
(A.29)
The lumped parameters of each cell are calculated by: hi hi _ in hi _ out 2 i f hi , P Ti f hi , P
(A.30)
(1) Cell 1 The density of working fluid at the outlet of Cell 1:
1_out in Mass conservation of working fluid:
(A.31)
Journal Pre-proof dL A L1 1 1 1_ out 1 min m1_ out dt t
(A.32)
Energy conservation of working fluid:
dh d dP dL1 A 1 1 h1 1 min hin m1_ out h1_ out i1 Di Tw1 T1 L1 1h1 1_ out h1_ out dt dt dt dt (A.33) Energy conservation of the tube wall: Cw w Aw
dTw1 i1 Di T1 Tw1 o1 Do Ths1 Tw1 dt
(A.34)
(2) Cell i (1
i_o ut i _ in in i
(A.35)
Mass conservation of working fluid: d Li i dL A Li i _ in 1_ out i 1 i i_ out i t dt dt
m m i _ in i_ out
(A.36)
Energy conservation of working fluid:
d i dhi i dP A i hi hi 1 Li dt dt P dt d Li dL A i _ in hi _ in i_ out hi_ out i 1 A i hi i _ out hi _ out i dt dt mi _ in hi _ in mi_ out hi_ out ii Di Twi Ti
(A.37)
Energy conservation of the tube wall: Cw w Aw
dTwi ii Di Ti Twi o1 Do Thsi Twi dt
(3) Cell N The heat transfer length of Cell N:
(A.38)
Journal Pre-proof LN L Li
(A.39)
N 1
Mass conservation of working fluid: d Li A N _ in N N 1 LN N dt P mN _ in mout
1 N 2 h h
P
dhN _ in dP 1 LN N dP dt 2 h
P
dhout dt
(A.40) Energy conservation of working fluid:
A N _ in hN _ in N hN
d Li N 1
dt
N 1 ALN N hN 2 hout
dhout dt P dP hN N 1 P h dt
1 dh dhN _ in 1 ALN N N _ in hN N 2 dP 2 h P dP mN _ in hN _ in mout hout iN Di LN TwN TN
(A.41)
Energy conservation of the tube wall: Cw w Aw
dTwN iN Di TN TwN oN Do ThsN TwN dt
(A.42)
Appendix B: Equations of MB model of the condenser. (1). Superheated region: Eqs.(A,1)-(A.3) are integrated over the subcooled region from 0 to L0 as shown in Fig.6 and are simplified with Leibniz’s rule. The lumped parameters of superheated region are calculated by:
hsuper hin hv 2 super f hsuper , P Tsuper f hsuper , P
(A.43)
Journal Pre-proof Mass conservation of working fluid:
dL 1 super A super v 0 L0 super P dt 2 hv h
P
dhv dP 1 super L0 dP dt 2 hin
P
dhin min m1 dt (A.44)
Energy conservation of working fluid: A super hsuper v hv
dL0 1 A L0 super hsuper super dt 2 hin
1 dh 1 AL0 super v hsuper super 2 dP 2 hv
P
dhv hsuper dP
min hin m1hv i 0 Di L0 Tw0 Tsuper
dhin dt P dP super 1 P h dt
(A.45)
Energy conservation of the wall:
dT T T dL Cw w Aw w0 w0 w1 0 i 0 Di T0 Tw0 o 0 Do Ths 0 Tw0 L0 dt dt
( A.46)
(2). Two-phase region: Eqs.(A,1)-(A.3) are integrated from L0 to L0 L1 . The lumped parameters of two-phase region are calculated by:
tp v 1 l tp htp v hv 1 l hl
(A.47)
Mass conversation of working fluid:
dL d dP dL d v A v l 0 v l 1 L1 1 l m1 m2 dt dt dP dt dP (A.48) Energy balance equation of working fluid: d l hl dP d v hv dL dL A l hl v hv 0 1 l hl v hv 1 L1 1 1 dt dt dP dP dt
m1hv m2 hl i1 Di L1 Tw1 T1
Journal Pre-proof (A.49) Energy balance equation of the wall: Cw w Aw
dTw1 i1 Di T1 Tw1 o1 Do Ths1 Tw1 dt
(A.50)
(3) Subcooled region: Eqs.(A,1)-(A.3) are integrated from L0 L1 to L. The lumped parameters of superheated region are calculated by:
hsub hout hl 2 sub f hsub , P Tsub f hsub , P
(A.51)
Mass balance equation of working fluid:
dL dL 1 sub A l sub 0 l sub 1 L2 sub P h 2 hl dt dt m2 mout
P
dhl dP 1 L2 dP dt 2 hl
P
dhl dt
(A.52) Energy balance equation of working fluid:
A l hl sub hsub
d L0 L1 1 AL2 sub hsub sub dt 2 hout
dhout dt P
1 dP dh 1 dhl AL2 sub l hsub sub hsub sub 1 dP 2 h '' P dP P h dt 2 m2 hl mout hout i 2 Di L2 Tw 2 T2
(A.53)
Energy balance equation of the wall:
dT T T dL dL Cw w Aw w 2 w1 w 2 1 2 i 2 Di T2 Tw 2 o 2 Do Ths 2 Tw 2 L2 dt dt dt (A.54)
Reference [1] Li J, Ge Z, Duan Y, Yang Z. Design and performance analyses for a novel
Journal Pre-proof organic Rankine cycle with supercritical-subcritical heat absorption process coupling. Applied Energy. 2019;235:1400-14. [2] Vetter C, Wiemer H-J, Kuhn D. Comparison of sub- and supercritical Organic Rankine Cycles for power generation from low-temperature/low-enthalpy geothermal wells, considering specific net power output and efficiency. Applied Thermal Engineering. 2013;51(1-2):871-9. [3] Yagli H, Koc Y, Koc A, Gorgulu A, Tandiroglu A. Parametric optimization and exergetic analysis comparison of subcritical and supercritical organic Rankine cycle (ORC) for biogas fuelled combined heat and power (CHP) engine exhaust gas waste heat. Energy. 2016;111:923-32. [4] Schuster A, Karellas S, Aumann R. Efficiency optimization potential in supercritical Organic Rankine Cycles. Energy. 2010;35(2):1033-9. [5] Dai X, Shi L, An Q, Qian W. Screening of hydrocarbons as supercritical ORCs working fluids by thermal stability. Energy Conversion and Management. 2016;126:632-7. [6] Chen Y, Lundqvist P, Johansson A, Platell P. A comparative study of the carbon dioxide transcritical power cycle compared with an organic rankine cycle with R123 as working fluid in waste heat recovery. Applied Thermal Engineering. 2006;26(17):2142-7. [7] Mikielewicz D, Mikielewicz J. A thermodynamic criterion for selection of working fluid for subcritical and supercritical domestic micro CHP. Applied Thermal Engineering. 2010;30(16):2357-62. [8] Xu H, Gao N, Zhu T. Investigation on the fluid selection and evaporation parametric optimization for sub- and supercritical organic Rankine cycle. Energy. 2016;96:59-68. [9] Chen H, Goswami DY, Rahman MM, Stefanakos EK. A supercritical Rankine cycle using zeotropic mixture working fluids for the conversion of low-grade heat into power. Energy. 2011;36(1):549-55. [10] Radulovic J, Castaneda NIB. On the potential of zeotropic mixtures in supercritical ORC powered by geothermal energy source. Energy Conversion and Management. 2014;88:365-71. [11] Dong B, Xu G, Cai Y, Li H. Analysis of zeotropic mixtures used in high-temperature Organic Rankine cycle. Energy Conversion and Management. 2014;84:253-60. [12] Xi H, Li M-J, He Y-L, Zhang Y-W. Economical evaluation and optimization of organic Rankine cycle with mixture working fluids using R245fa as flame retardant. Applied Thermal Engineering. 2017;113:1056-70. [13] Xie H, Yang C. Dynamic behavior of Rankine cycle system for waste heat recovery of heavy duty diesel engines under driving cycle. Applied Energy. 2013;112:130-41. [14] Quoilin S, Aumann R, Grill A, Schuster A, Lemort V, Spliethoff H. Dynamic modeling and optimal control strategy of waste heat recovery Organic
Journal Pre-proof Rankine Cycles. Applied Energy. 2011;88(6):2183-90. [15] Horst TA, Tegethoff W, Eilts P, Koehler J. Prediction of dynamic Rankine Cycle waste heat recovery performance and fuel saving potential in passenger car applications considering interactions with vehicles’ energy management. Energy Conversion and Management. 2014;78:438-51. [16] Zhang J, Zhang W, Hou G, Fang F. Dynamic modeling and multivariable control of organic Rankine cycles in waste heat utilizing processes. Computers & Mathematics with Applications. 2012;64(5):908-21. [17] 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 Conversion and Management. 2017;148:724-36. [18] Chowdhury JI, Bao Kha N, Thornhill D. Dynamic model of supercritical organic rankine cycle waste heat recovery system for internal combustion engine. International Journal of Automotive Technology. 2017;18(4):589-601. [19] Chowdhury JI, Bao Kha N, Thornhill D, Hu Y, Soulatiantork P, Balta-Ozkan N, et al. Fuzzy Nonlinear Dynamic Evaporator Model in Supercritical Organic Rankine Cycle Waste Heat Recovery Systems. Energies. 2018;11(4). [20] Wang T, Gao N, Zhu T. Investigation on the optimal condensation temperature of supercritical organic Rankine cycle systems considering meteorological parameters. Energy Conversion and Management. 2018;174:54-64. [21] Olumayegun O, Wang M. Dynamic modelling and control of supercritical CO2 power cycle using waste heat from industrial processes. Fuel. 2019;249:89-102. [22] He C, Liu C, Gao H, Xie H, Li Y, Wu S, et al. The optimal evaporation temperature and working fluids for subcritical organic Rankine cycle. Energy. 2012;38(1):136-43. [23] Jian L, Qiang L, Zhong G, Duan Y, Zhen Y. Thermodynamic performance analyses and optimization of subcritical and transcritical organic Rankine cycles using R1234ze(E) for 100–200 °C heat sources. Energy Conversion & Management. 2017;149:140-54. [24] Lemmon EW, Huber ML, Mclinden MO. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 9.1. NIST NSRDS -. 2010. [25] Xu XH, Duan YY, Yang Z. Prediction of the critical properties of binary alkanol + alkane mixtures using a crossover CPA equation of state. Fluid Phase Equilibria. 2011;309(2):168-73. [26] Yang F, Liu Q, Duan Y, Yang Z. Crossover multiparameter equation of state: General procedure and demonstration with carbon dioxide. Fluid Phase Equilibria. 2019;494:161-71. [27] Kiselev SB. Cubic crossover equation of state1Contribution of the National Institute of Standards and Technology, not subject to copyright in the United States, November 26, 1997.1. Fluid Phase Equilibria. 1998;147(1):7-23. [28] Xu XH, Duan YY. Crossover CPA equation of state for associating fluids.
Journal Pre-proof Fluid Phase Equilibria. 2010;290(1):148-52. [29] Dicko M, Coquelet C. Application of a new crossover treatment to a generalized cubic equation of state. Fluid Phase Equilibria. 2011;302(1):241-8. [30] Shen A-J, Liu Q, Duan Y-Y, Yang Z. Crossover VTSRK equation of state for selected alkane + alkane and CO2 + alkane binary mixtures. Fluid Phase Equilibria. 2016;408:180-9. [31] Belyakov MY, Kulikov VD, Muratov AR, Voronov VP. Crossover Equation of State of a Multi-Component Fluid Mixture in The Vicinity of Liquid-Vapor Critical Points. Chemical Physics. 2018;513:S0301010418303173. [32] Oyewunmi OA, Taleb AI, Haslam AJ, Markides C. An assessment of working-fluid mixtures using SAFT-VR Mie for use in organic rankine cycle systems for waste-heat recovery. 2014;6(4):301-16. [33] Gnielinski V. New equations for heat and mass transfer in turbulent pipe and channel flow. International Chemical Engineering. 1976;16:359-68. [34] Deutsches MO, Elmqvist H, Ab D, Otter M, Jaschinski A, Schweiger C, et al. Moving Boundary Models for Dynamic Simulations of Two-phase Flows. The second international modelica conference. Germany, March 18-19; 2002. [35] Zivi SM. Estimation of Steady-State Steam Void-Fraction by Means of the Principle of Minimum Entropy Production. Transasme Serc. 1964;86(2):247. [36] Bell KJ, Ghaly MA. An Approximate Generalized Design Method for Multicomponent/Partial Condensers. AIChE Symp Ser. 1973;69:72-9. [37] Shah MM. An Improved and Extended General Correlation for Heat Transfer During Condensation in Plain Tubes. HVAC&R Research. 2009;15(5):889-913. [38] Le VL, Kheiri A, Feidt M, Pelloux-Prayer S. Thermodynamic and economic optimizations of a waste heat to power plant driven by a subcritical ORC (Organic Rankine Cycle) using pure or zeotropic working fluid. Energy. 2014;78:622-38. [39] Zhang Y, Deng S, Zhao L, Lin S, Bai M, Wang W, et al. Dynamic test and verification of model-guided ORC system. Energy Conversion and Management. 2019;186:349-67. [40] Jiménez-Arreola M, Pili R, Wieland C, Romagnoli A. Analysis and comparison of dynamic behavior of heat exchangers for direct evaporation in ORC waste heat recovery applications from fluctuating sources. Applied Energy. 2018;216:724-40. [41] Dalkilic AS, Yildiz S, Wongwises S. Experimental investigation of convective heat transfer coefficient during downward laminar flow condensation of R134a in a vertical smooth tube. International Journal of Heat and Mass Transfer. 2009;52(1):142-50.
Journal Pre-proof Highlights
The improved dynamic model with higher accuracy and applicability is developed. The effects around pseudo-critical point of working fluid are investigated. A prediction function of response time and heat transfer coefficients is found. The effects of heat source frequencies and system thermal inertia are researched.