Fuzzy Nonlinear Dynamic Evaporator Model in Supercritical Organic Rankine Cycle Waste Heat Recovery Systems

Abstract

The organic Rankine cycle (ORC)-based waste heat recovery (WHR) system operating under a supercritical condition has a higher potential of thermal efficiency and work output than a traditional subcritical cycle. However, the operation of supercritical cycles is more challenging due to the high pressure in the system and transient behavior of waste heat sources from industrial and automotive engines that affect the performance of the system and the evaporator, which is the most crucial component of the ORC. To take the transient behavior into account, the dynamic model of the evaporator using renowned finite volume (FV) technique is developed in this paper. Although the FV model can capture the transient effects accurately, the model has a limitation for real-time control applications due to its time-intensive computation. To capture the transient effects and reduce the simulation time, a novel fuzzy-based nonlinear dynamic evaporator model is also developed and presented in this paper. The results show that the fuzzy-based model was able to capture the transient effects at a data fitness of over 90%, while it has potential to complete the simulation 700 times faster than the FV model. By integrating with other subcomponent models of the system, such as pump, expander, and condenser, the predicted system output and pressure have a mean average percentage error of 3.11% and 0.001%, respectively. These results suggest that the developed fuzzy-based evaporator and the overall ORC-WHR system can be used for transient simulations and to develop control strategies for real-time applications.

1. Introduction

Around 60% of global greenhouse gases are emitted to the environment from industry and transport sectors [1]. These two sectors are heavily involved in energy conversion processes such as the burning of fuels to generate electricity, heat, and mechanical power. Low efficiency in the conversion process is one of the determinants that causes increased pollutant emissions. In a typical diesel internal combustion engine, a maximum of 45% of fuel energy can be converted into the mechanical energy at its best operating condition, while the gasoline engine returns a high of 35% [2]. Remaining energy is wasted mainly through heat lost to the engine’s coolant and exhaust. To mitigate the environmental impact caused by the higher fuel consumption and the lower conversion efficiency, recent research has been focused on the development and utilization of alternative sources of energy such as renewable energies, different types of waste heat and biomass. Waste heat from industries and automotive sectors has gained much attention to improve the energy conversion efficiency of the engine and reduce the environmental pollutions. Since the waste heat from these sectors are mostly low to medium grade in quality, one way of utilizing this heat is to recover and convert into mechanical rotations or electrical power using an appropriate waste heat recovery (WHR) system, which increases the engine’s thermal efficiency and thus reduces the fuel consumptions and emissions.

Among different WHR technologies, thermoelectric generator (TEG), phase change material (PCM) engine, Stirling engine, and organic Rankine cycle (ORC) are most common. A TEG can convert the heat energy to electrical power when a temperature difference is applied to a joint of two different materials. A PCM engine generates mechanical rotation using the expansion and contraction properties of PCMs used for WHR. A Stirling engine produces mechanical power through a continuous change of volume of a gas used inside the cylinders of the engine. The waste heat in an ORC is used to vaporize an organic working fluid, elevating its pressure and temperature and then expanding the fluid to generate mechanical energy. The TEG, PCM and Stirling engines are less flexible to construct a system for low to medium grade WHR with a kW range power output; while an ORC has a great flexibility to choose different components suitable for many ranges of power output. Moreover, the conversion efficiency of the PCM, TEG, Stirling Engine, and the ORC in low to medium grade WHR applications are, generally, up to 2.5%, 5%, 10.3%, and 16%, respectively [3,4]. Although some reports suggest that Stirling engines have many advantages such as low noise and high thermal efficiency for medium to high grade heat recovery applications, this technology is still considered to be in an early stage of development [5]. On the contrary, the ORC has many practical advantages including adapting with different grades of heat sources and availability of its components as reported in [3,6,7].

The operating cycle of an ORC can be one of the two types: subcritical cycle whose pressure is lower than the critical pressure of the organic fluid or supercritical cycle whose pressure is higher than the critical pressure of the fluid (see Figure 1). The cycle $1–2–2^{\prime }–3–4–5–5^{\prime }–1$ represents a subcritical ORC, while $1–2–3^{\prime }–5^{\prime }–1$ represents a supercritical ORC. The subcritical ORC for WHR applications is widely used and mentioned in the following reports [8,9,10,11,12,13,14]. Authors in [15] concluded that in subcritical conditions the thermal efficiency of the cycle is low because of higher exergy losses and destructions while the authors in [6,7,15,16] showed that these losses and destructions become lower in supercritical conditions that lead to high thermal efficiency of the ORC. Moreover, the cycle efficiency of the supercritical ORC is also dependent on system configurations, types of working fluids and heat sources used in the system. Lecompt et al. [17] compared the efficiency of low temperature subcritical and supercritical WHR cycles and reported that the supercritical cycle outperforms the subcritical by 10.8% in cycle efficiency. A similar investigation carried out by Chen et al. [18] indicated that up to 30% increase in the cycle efficiency is achievable if a WHR system is run at a supercritical pressure rather than a traditional subcritical pressure.

The analysis and performance evaluation of supercritical ORC cycles mainly focused on fluid selections [6,15,19,20], design and optimization [15,21,22,23,24,25,26,27]. Although steady-state models of ORC-WHR systems are necessary for the preliminary analyses as discussed above, they cannot be used either for performance evaluation in transient conditions or control system simulations. Therefore the dynamic model of the ORC-WHR system is proposed by the authors in [28]. Since the dynamics of the overall ORC systems evolve around the evaporator, the critical component of the system, the model is developed using conventional finite volume (FV) method. In this method, the evaporator is discretized into several finite segments and heat transfer equations for each segment are solved numerically. The FV is a robust and accurate method which is commonly used as can be seen in [12,29,30]. Despite the fact that the FV model is accurate and able to capture the variable thermo-physical properties of the supercritical fluids as described in [31], it is very time intensive for computation and this limits its real-time applications. To eliminate the limitation, a novel fuzzy-based dynamic evaporator mode is developed and presented in this paper. The knowledge and transient response of the evaporator and the ORC system under variable heat input conditions analyzed using the FV model is used to predict the thermal inertia of the fuzzy model. Other components of the cycle are also developed and then integrated with the novel fuzzy-based dynamic evaporator model to investigate the performance of the system under a variable heat source. The detailed performance of the evaporator model, which is based on fuzzy and the complete ORC-based WHR system are validated and presented in this research.

This paper is organized as follows; the next section presents the system description of the ORC-WHR considered in this research. Section 3 presents the individual components models including FV and fuzzy-based evaporator, model validations and integrations. The performance of the dynamic system is evaluated and discussed in Section 4. Finally, the research findings and conclusions are presented in Section 5.

2. ORC-WHR System

The ORC-WHR system studied consists of a pump, an evaporator, an expander couple with a generator, a condenser, and an accumulator, as shown in Figure 2. Additionally, a valve in dynamic conditions is used to manage the system start-up and shut-down procedures. The ORC is similar to the conventional steam Rankine cycle but uses organic fluids instead of water as a working fluid. The organic fluid has desirable thermo-physical properties, such as low boiling point and high vapor pressure, that are suitable for low to medium grade heat recovery applications with a output capacity from a few kilowatts to a couple of hundred kilowatts [8,29,32,33]. In this research, R134a was chosen as the working fluid as it has advantages such as availability, wide commercial use, low critical temperature, high auto-ignition temperature, and suitability for low grade heat recovery applications. The thermo-physical properties of R134a are shown in

Table 1. Properties of R134a.

Properties	Value
Chemical formula	CH2FCF3
Molar mass	102.03 kg/kmol
Critical point density	511.9 kg/m3
Critical Pressure	4.06 MPa
Maximum pressure	70 MPa
Boiling point temperature	247.08 K
Critical temperature	374 K
Autoignition temperature	1016 K

. The working fluid of the ORC is pumped to the evaporator and heated up by a heat source to its evaporating point, e.g., a pressurized hot water, and then expanded in the expander to drive a coupled generator. The excess heat contained in the working fluid is released at the condenser, and then the cold fluid is sent back to the accumulator and completes the cycle. The ORC-WHR system will be simulated at a pressure of 6 MPa, which is a supercritical pressure that is much higher than the working fluid used in most thermal cycles.

3. Model Development

The components of the ORC-WHR system described in the previous section are modelled and integrated in this section. The modelling of the valve, however, is disregarded due to that fact that the start-up and shut-down procedures are not covered in this paper.

3.1. Evaporator Model

The evaporator is a critical component of the ORC system since effective heat transfer from this component influences the thermal efficiency of the system. Also, due to the transient nature of heat sources from industry and automotive sectors, the evaporator is primarily treated as a source of thermal inertia in the system. Several heat exchangers can be employed as the evaporator in WHR systems. These are the finned tube, the shell and tube and the plate heat exchangers, etc. Plate heat exchangers are highly recommended in this application because they are compact and have a large area, which aims to recover the most possible quantity of heat from the main heat source [5]. Other advantages of plate heat exchangers are also a minimum risk of internal leakage, minimum pressure drop, and less complexity in maintenance [34].

3.1.1. Finite Volume Model

The FV method is used to develop the dynamic model of the evaporator where the heat exchanger is discretized into a finite number of segments along the direction of the flow (See Figure 3). The number of segments, N on both sides and wall of the heat exchanger, is the same. In the numerical model, the input and output parameters of the working fluid and the heat source are imposed on the inlet and outlet sides of each segment as shown in Figure 3. These parameters are stored at locations called nodes. The number of nodes n is therefore, n = N + 1.

Below assumptions are used to develop the FV evaporator model:

• The heat exchanger model is assumed to be one-dimensional, and the heat transfer to surrounding environment is neglected.
• The momentum conservation is not considered in the model and the pressure variation within the model is assumed to be negligible.
• The heat transfer between the heat source and the refrigerant takes place not by conduction but by convection.
• Heat exchanger wall is uniformly built, and thermo-physical properties are assumed to be constant.
• Thermo-physical properties of refrigerant and heat source fluid for each discrete segment are constant.

All these assumptions are commonly used to simplify the modelling effort and reduce the computational time, which can be seen in the literature [12,22,29,35]. Details of mass and energy conservation equations of the FV method are completely discussed by the authors in [28].

Solution Methodology

In the FV approach, the heat source and the refrigerant outlet conditions have been initially guessed. Using the initial and inlet conditions, the simultaneous first-order differential equations for the wall, the refrigerant, and heat source of the FV model described in [28] are solved accordingly. The time-dependent terms in the governing equations are integrated using the Euler explicit method. This method is one of the known approaches used for the solution of time-dependent ordinary or partial differential equations. In this approach, the numerical approximation of the time-dependent term at the next time step is calculated from the state of the term at the current step. The numerical discretization of the FV evaporator model in the time domain can be represented as follows:

$$H_{r,j}{}^{t+1}=H_{r,j}{}^t+\frac{\Delta t}{V_{r,j}{\rho }_{r,j}{}^t}\left(A_{r,j}{h}_{r,j}{}^t(T_{wall,j}{}^t-T_{r,j}{}^t)-{\dot{m}}_{r,j}{}^t(H_{r,j}{}^t-H_{r,j-1}{}^t)\right)$$

(1)

$$T_{h,j}{}^{t+1}=T_{h,j}{}^t+\frac{\Delta t}{V_{h,j}{\rho }_{h,j}{}^t{c}_p{{}_{h,j}}^t}\left({\dot{m}}_{h,j}{}^t{c}_p{{}_{h,j}}^t(T_{h,j+1}{}^t-T_{h,j}{}^t)-A_{h,j}{h}_{h,j}{}^t(T_{h,j}{}^t-T_{wall,j}{}^{t+1})\right)$$

(2)

$$T_{wall,j}{}^{t+1}=T_{wall,j}{}^t+\frac{\Delta t}{m_{wall,j}{c}_p{}_{wall,j}}\left(A_{r,j}{h}_{h,j}{}^t(T_{h,j}{}^t-T_{wall,j}{}^t)-A_{r,j}{h}_{r,j}{}^t(T_{wall,j}{}^t-T_{r,j}{}^t)\right)$$

(3)

where $h_h$ and $h_r$ are the convective heat transfer coefficients (kW/m²K) of the heat source fluid and refrigerant with the wall; $c_p{}_h$ and $c_p{}_{wall}$ (kJ/kgK) are the specific heat capacity of the heat source and heat exchanger wall, respectively. ${\dot{m}}_r$ (kg/s), $T_r$ (K) and $H_r$ (KJ/kg K) are the mass flow rate, temperature and enthalpy of the refrigerant; ${\dot{m}}_h$, $T_h$ and $T_{wall}$ are the mass flow rate, temperature of heat source and temperature of wall, respectively.

Time-Step Determination

In order to achieve a stable solution of the discretized Equations (1)–(3), a correlation to calculate the time step $\Delta t$ can be derived from the equations that contain the variables of interest $H_{r,j}{}^t$ and $T_{h,j}{}^t$. To obtain a stable solution, the coefficients of $H_{r,j}{}^t$ and $T_{h,j}{}^t$ variables in Equations (1)–(3) should be greater than zero. Therefore, from Equation (1), it can be written as follows:

$$\Delta t\le \frac{V_{r,j}{\rho }_{r,j}{}^t}{{\dot{m}}_{r,j}{}^t}$$

(4)

Similarly, from Equation (2), this can be written as,

$$\Delta t\le \frac{V_{h,j}{\rho }_{h,j}{}^t{c}_{p_{h,j}}{}^t}{{\dot{m}}_{h,j}{}^t{c}_{p_{h,j}}{}^t+A_h{h}_{h,j}{}^t}$$

(5)

The evaporator model time step should be the minimum of the two values in Equations (4) and (5) which can be represented by Equation (6). The minimum time step in the discretized equations can guarantee a stable solution in the numerical model.

$$\Delta t=\min \left[\frac{V_{r,j}{\rho }_{r,j}{}^t}{{\dot{m}}_{r,j}{}^t},\frac{V_{h,j}{\rho }_{h,j}{}^t{c}_{p_{h,j}}{}^t}{{\dot{m}}_{h,j}{}^t{c}_{p_{h,j}}{}^t+A_h{h}_{h,j}{}^t}\right]$$

(6)

Simulation Procedure

The dynamic evaporator model consists of two parts as shown in Figure 4. They are the steady state part and the dynamic part. In the steady state calculation, the model is converged to the steady state condition before applying any transient inputs to the model. This convergence reduces the effect of initial assumptions to the accuracy when the model is simulated in the transient conditions. The steps involved in the solution of discretized Equations (1)–(3) are as follows:

Step 1: Define the geometrical parameters, heat source fluid and refrigerant inputs and a discrete time step of the evaporator including the number of segments and nodes of the model. REFPROP [36] database has been utilized to obtain the thermo-physical properties of the refrigerant and the heat source fluid during the simulation.

Step 2: The second phase is the initialization which is an important issue for every dynamic system in numerical methods. Given the steady inputs, the unknown parameters of Equations (1)–(3) are initialized by assuming the appropriate values at $t=0$ as follows:

$$H_{r,j}{}^t=H_{r,ini}{}^{t=0},T_{h,j}{}^t=T_{h,ini}{}^{t=0}$$

(7)

From the refrigerant’s initial enthalpy distribution, the temperature of the segment can be retrieved as follows:

$$T_{r,j,ini}{}^{t=0}=f(P_{ev},H_{r,ini}{}^{t=0})$$

(8)

The initial temperatures of the refrigerant and the heat source fluid are used to define the initial temperature of the wall as follows:

$$T_{wall,j}{}^t=T_{wall,j,ini}{}^{t=0}=\frac{T_{r,j,ini}{}^{t=0}+T_{h,j,ini}{}^{t=0}}{2}$$

(9)

The initial temperatures of the fluids and wall mass are assigned equally to all segments to have an initial temperature distribution for the model. This initial distribution does not necessarily have to be accurate since they serve only as a starting point for the simulation. Once all the parameters of the current time $t=0$ are available, the solution of the discretized Equations (1)–(3) can be started.

Step 3: Based on the enthalpy, temperature distribution and the initial inputs condition at the previous step, this step retrieves thermodynamic properties from REFPROP and calculates the convective heat transfer coefficients of the heat source fluid $h_{h,j}{}^t$ and refrigerant $h_{r,j}{}^t$ at $t=0$, respectively. Thermo-physical properties and the local temperature are considered in calculating the convective heat transfer coefficients in Equations (1)–(3), using the Nusselt number ($Nu$) and Reynolds number ($\mathrm{Re}$) correlations in Equation (10). The specific Nusselt number correlations for the heat source fluid and the refrigerant are discussed by the authors and can be found in [28].

$$Nu=\frac{hD}{K},\mathrm{Re}=\frac{\rho vD}{\mu }$$

(10)

where $D$, $\rho$, $\mu$ and $v$ are the hydraulic diameter of the plate heat exchanger, the density of the fluid in (kg/m³) and the viscosity in (Pa.s) and the velocity of the fluids (m/s), respectively.

Within this simulation step, the enthalpy, heat source temperature and wall temperature for the next time step are calculated using Equations (1)–(3). This simulation step is repeated until the calculation reaches to the final segment of the model.

Step 4: At this step of the simulation, the Nth values of the refrigerant enthalpy, heat source temperature and wall temperature at $t=t+1$ and at $t=0$ are compared. If the differences between consecutive values are less than the predefined steady-state convergence values of ${\delta }_{1,}{\delta }_2$ and ${\delta }_3$, which is 0.0001, the solution is said to have achieved the steady state condition. Otherwise, steps 1–3 are repeated until the steady state condition is reached. At the end of Step 4, the initial distribution of the refrigerant enthalpy, the heat source temperature, and the wall temperature are all converged to the actual value and the effect of the initial assumptions are now completely vanished.

Step 5: The dynamic part of the model is started at the end of the steady state part as depicted in Figure 4. The transient inputs to the model are introduced at this step. At a discrete time-step $\Delta t$ and model time $t=\Delta t_1$, the calculation of the state variables is continued along the segments $1<j<N$. At the end of each time step, the values of the state variables are stored in the MATLAB workspace. The dynamic simulation is continued until it reaches the final predefined model time $t=t_n$.

Numerical Issues

The discrete time step in Equation (6) depends on the mass flow rate of each fluid, the density of each fluid and the volume of each discrete segment. Since the mass flow rate and density of the fluids are time dependent variables and can vary during the simulation, the discrete time step should also vary according to Equation (6). Calculating a variable time step during the simulation will make the computation process more complicated and time-consuming. Moreover, this process has a limitation in a real-time application as the temperature at the evaporator inlet, and therefore the properties of the fluid is unknown and cannot be predicted in advance. Therefore, instead of defining a variable time step, a fixed time step is used based on the maximum and minimum temperature and flow rate ranges used in the simulation. In this process, the discrete time step is set to 0.1 s in the simulation, which is selected using Equation (6) that can guarantee a stable solution of the discretized Equations (1)–(3).

In the FV technique, the number of finite cells is also an important issue that influences the accuracy as well as the time requirement of the simulation. To evaluate the error of the calculation, the outputs of the model at different numbers of segments (10, 20 and 50) are compared with that of a 100 segments model, which is taken as the reference value. The computation errors correspond to 10, 20, and 50 segments model for the prediction of the refrigerant enthalpy are approximately 7%, 3.1%, and 0.8%, respectively. On the contrary, the calculation error for the outlet temperature of the heat source is below 1% for all number of segments in the simulation. It is as expected that the higher accuracy is achieved with the higher number of segments. However, this higher number of segments leads to an extensive computation time. Therefore, in this simulation the number of segment is set to 20, as it is a good trade-off between the model accuracy and the simulation time.

Transient Response Simulation

To evaluate the performance of the evaporator model in transient situations, three different tests have been carried out. The step change to the inlet temperature and mass flow rate of the heat source are the first and second test. Finally, the step change in the mass flow rate of the refrigerant is also investigated in this section.

Figure 5 depicts the evaporator model response to a step change in the heat source inlet temperature. The step magnitude of 48 K at t = 150 s, which is arbitrarily chosen, is applied to the inlet temperature of the heat source. To discard the effects of the assumptions initially considered, the model is converged into the steady state condition before the step change occurred at t = 150 s. The other variables at the inlet are kept constant during this transient response test, and the calculation is continued until the steady state condition is obtained.

To investigate the effect of a mass flow rate transient on the model outputs, a step change of 0.2 kg/s to the mass flow rate of the heat source is imposed at t = 150 s. Figure 6 shows the responses of refrigerant enthalpy, heat source temperature and wall temperature to the step change. Similarly, a step change of 0.15 kg/s to the refrigerant flow rate is also imposed and the transient responses of the model are illustrated in Figure 7. The outlet variable responses due to the step changes in the heat source and refrigerant inputs indicate that the proposed model is stable and can be used to simulate the transient behavior of the inlet conditions.

The response time of the evaporator model to different transient steps is shown in

Table 2. Ideal response time of finite volume (FV) evaporator model outputs.

Constant Input Conditions	Transient Condition (Step Change)	Ideal 63% Response Time	Enthalpy and Temperature at 63% Response Time
${\dot{m}}_r=0.1$ kg/s	$T_h=48$ K Step down	${\zeta }_{H_{\_r\_o}}=16.6s$	$H_{\_r\_o,}{}_{\zeta }=431.94$ kJ/kg
$T_r=303$ K
${\dot{m}}_h=0.3$ kg/s		${\zeta }_{T_{\_h\_o}}=27.1s$	$T_{\_h\_o,}{}_{\zeta }=478.96$ K
${\dot{m}}_r=0.1$ kg/s	${\dot{m}}_h=0.2$ kg/s Step down	${\zeta }_{H_{\_r\_o}}=23.7s$	$H_{\_r\_o,}{}_{\zeta }=447.1$ kJ/kg
$T_r=303$ K
$T_h=523$ K		${\zeta }_{T_{\_h\_o}}=53.5s$	$T_{\_h\_o,}{}_{\zeta }=490.5$ K
${\dot{m}}_h=0.3$ kg/s	${\dot{m}}_r=0.15$ kg/s Step down	${\zeta }_{H_{\_r\_o}}=25.5s$	$H_{\_r\_o,}{}_{\zeta }=436.6$ kJ/kg
$T_r=303$ K
$T_h=523$ K		${\zeta }_{T_{\_h\_o}}=16.8s$	$T_{\_h\_o,}{}_{\zeta }=502.98$ K

. The response time is a time that is required for a time-dependent variable to reach from one steady state value to 63% of the final steady state value [37]. The state values of the evaporator model at the 63% response time are also shown in Table 2. It can be seen from Figure 5, Figure 6 and Figure 7 and Table 2 that the ideal response time is different with respect to different step changes, which is as expected since the modification in one inlet condition will change all outlet conditions according to their nonlinear relationship as represented by Equations (1)–(3). The results from the transient response analysis indicate that the behavior of the evaporator model is satisfactory and therefore can be used in the simulation of the system dynamics.

Simulation Time Constraint in Dynamic Scenario

There are several factors that can increase the simulation time of the evaporator model in both the steady state and dynamic conditions. The first factor is the discrete time step and the second factor contributes to the higher simulation time is the number of segments of the model.

In the dynamic situation, since the model initially needs to converge into a steady state condition before applying the transient heat source, the simulation time of the model can be evaluated by combining the time required for both conditions. The simulation time of the evaporator model with the random transient heat source in Figure 8 is presented in

Table 3. Computation time for a transient heat source simulation.

No. of Segments	Total Simulated Run-Time of the System (Steady + Dynamic)	Actual Computing Time of the Simulation
10	1300 s	2040 s
20	1300 s	3820.6 s
50	1300 s	9131.8 s
100	1300 s	17,782 s

. It can be concluded from Table 3 that despite the total simulated run-time of the system being kept the same, the actual computing time of the simulation is increased with the number of segments. The model with this higher simulation time will show poor performance in real-time control applications. Therefore, to reduce the simulation time, a second approach to develop the evaporator model using the fuzzy technique is proposed in the next section.

3.1.2. Fuzzy Based Dynamic Evaporator Model

Almost all dynamic systems in practical applications are nonlinear. Modelling of those systems cannot be represented by simple linear differential equations. It requires a set of high order nonlinear equations which make the whole system too complex to be used in a control system. In the WHR system, the evaporator operating under variable heat source is a nonlinear. Nonlinear models of the evaporator in a WHR system have been developed and are shown in several reports [30,38]. However, even with the simplified assumption and model order reduction made for the evaporator, high nonlinearity still makes the system complex and time-consuming for the simulation [39]. The simplified 1D FV model of the evaporator developed in this research is also a nonlinear model and required extensive time in simulations as shown in Table 3. The nonlinearity of the model comes from the state variables, which are the functions of the thermodynamic properties and input conditions according to Equations (1)–(3). For this reason, this section presents a dynamic model of the evaporator based on fuzzy logic. The fuzzy logic is a multivalued logical system that provides the value of an unknown output by attaching the degree of known input and output of the system [40]. The fuzzy logic concept can be used to develop a model of a nonlinear system which does not require a set of complex mathematical equations. Therefore, it requires less resource and saves substantial computation time.

The design of the dynamic fuzzy-based evaporator model is similar to that of the FV model. The evaporator model has three inputs (refrigerant mass flow rate ${\dot{m}}_r$, heat source mass flow rate ${\dot{m}}_h$ and heat source temperature $T_h$) and two outputs (outlet temperature of the refrigerant $T_{r,o}$ and outlet temperature of heat source $T_{h,o}$).

The dynamic fuzzy model is composed of two parts: the fuzzy inference system and the thermal inertia of the system as shown in Figure 9. The outputs of the fuzzy inference system, $T_{r,o}$ and $T_{h,o}$, are fed to the thermal inertia calculation block. The dynamic model outputs, $T_{\_r\_o}$ and $T_{\_h\_o}$, are then obtained from the thermal inertia block.

Based on the data available from the FV dynamic model, the range of the inputs and output parameters are determined and adjusted from the experience about the system of the designer as follows: $[{\dot{m}}_{r,\min },{\dot{m}}_{r,\max }]=$ [0.01 (kg/s), 0.25 (kg/s)], $[{\dot{m}}_{h,\min },{\dot{m}}_{h,\max }]=$ [0.05 (kg/s), 0.300 (kg/s)], $[T_{h,\min },T_{h,\max }]=$ [400 (K), 523 (K)], $[T_{r,o,\min },T_{r,o,\max }]=$ [364.5 (K), 470 (K)] and $[T_{h,o,\min },T_{h,o,\max }]=$ [285 (K), 595 (K)], respectively. To obtain a high efficiency of the fuzzy inference system, the input and output variables are normalized in between 0 and 1 as illustrated in Figure 10 and Figure 11. The actual crisp values of the model parameters are calculated from the normalized values as follows:

$${\dot{m}}_r=0.01+0.24{\dot{m}}_r{}^{\prime }$$

(11)

$${\dot{m}}_h=0.05+0.25{\dot{m}}_h{}^{\prime }$$

(12)

$$T_h=400+123T_h{}^{\prime }$$

(13)

$$T_{r,o}=364.5+105.5T_{r,o}{}^{\prime }$$

(14)

$$T_{h,o}=285+310T_{h,o}{}^{\prime }$$

(15)

where ${\dot{m}}_r{}^{\prime }$, ${\dot{m}}_h{}^{\prime }$, $T_h{}^{\prime }$, $T_{r,o}{}^{\prime }$ and $T_{h,o}{}^{\prime }$ are the normalized values of their corresponding variables.

The fuzzy-based model input and output ranges must be classified into different linguistic levels, which are called membership functions. The linguistic level assigned to the input variables ${\dot{m}}_r$, ${\dot{m}}_h$, $T_h$ are: L: Low; M: Medium and H: High. The output variable $T_{r,o}$ is designed by five membership functions: VL: Very Low; L: Low; M: Medium; H: High and VH: Very High; and $T_{h,o}$ has one additional membership function than $T_{r,o}$ that is Medium to High, MH. The numbers of membership functions were selected in a way so that they can capture the transient nature of the evaporator at an accuracy over 90%. All the inputs-outputs variables and their membership functions are shown in Figure 10 and Figure 11.

Once the membership functions of the input and output variables are all set, the next step is to define the rules and create a logical relationship among the membership functions of the model inputs and the model outputs. The fuzzy rules are established using the abovementioned sets of the input and output variables as follows:

$$\mathrm{Rule}i:\mathrm{IF}{\dot{m}}_r\mathrm{is}{\alpha }_i\mathrm{AND}{\dot{m}}_h\mathrm{is}{\beta }_i\mathrm{AND}{T}_h\mathrm{is}{\gamma }_i\mathrm{THEN}{T}_{r,o}\mathrm{is}{\delta }_i\mathrm{AND}{T}_{h,o}\mathrm{is}{\psi }_i$$

(16)

where $i=1,2,3.....n$, $n$ is the number of fuzzy rules, ${\alpha }_i$, ${\beta }_i$, ${\gamma }_i$, ${\delta }_i$, ${\psi }_i$ are the $i$th fuzzy sets of the input and output variables of the fuzzy system.

The rules of the fuzzy model are generally defined from the knowledge of characteristics of a system. In this dynamic model, rules are established from the intuition and knowledge of the evaporator used in this paper. The rules of the fuzzy dynamic model are presented in

Table 4. Fuzzy rules of dynamic model of the evaporator. VL: Very Low; L: Low; M: Medium; MH: Medium to High; H: High and VH: Very High.

Rule Number	IF ${\dot{\mathit{m}}}_{\mathit{r}}$ is	AND ${\dot{\mathit{m}}}_{\mathit{h}}$ is	AND ${\mathit{T}}_{\mathit{h}}$ is	THEN ${\mathit{T}}_{\mathit{r},\mathit{o}}$ is	AND ${\mathit{T}}_{\mathit{h},\mathit{o}}$ is
1	L	L	L	L	VL
2	L	L	M	M	L
3	L	L	H	M	L
4	L	M	L	L	L
5	L	M	M	L	LM
6	L	M	H	VH	M
7	L	H	L	M	LM
8	L	H	M	H	M
9	L	H	H	VH	MH
10	M	L	L	VL	VL
11	M	L	M	L	L
12	M	L	H	L	L
13	M	M	L	VL	LM
14	M	M	M	L	M
15	M	M	H	M	MH
16	M	H	L	VL	M
17	M	H	M	L	MH
18	M	H	H	M	VH
19	H	L	L	VL	VL
20	H	L	M	VL	L
21	H	L	H	VL	L
22	H	M	L	VL	LM
23	H	M	M	L	M
24	H	M	H	L	MH
25	H	H	L	VL	M
26	H	H	M	L	H
27	H	H	H	L	VH

. Overall mapping of the input parameters to the outputs represented by the fuzzy rules can be plotted as 3D surfaces. The fuzzy surfaces of the evaporator and heat source outlet temperature are shown in Figure 12 and Figure 13.

Among different fuzzy reasoning methods, the MAX-MIN method [41] is utilized to obtain the output from the present input and the inference rules. For a given particular input fuzzy set ${\Omega }^{\prime }$ in $U$, the output fuzzy set ${\Phi }^{\prime }$ in $S$ for $T_{r,o}$ is computed through the inference system is as follows:

$${\mu }_{\Phi }^{\prime }(T_{r,o})=\underset{l=1}{\overset{n}{\max }}\left[\underset{x\in U}{\sup }\min \left\{{\mu }_{{\Omega }^{\prime }}(x),{\mu }_{{\Omega }_1^l}({\dot{m}}_r),{\mu }_{{\Omega }_2^l}({\dot{m}}_h),{\mu }_{{\Omega }_3^l}(T_h),{\mu }_{{\Phi }^l}(T_{r,o})\right\}\right]$$

(17)

The output membership functions for $T_{h,o}$ is calculated similarly. In this research, the membership functions used for the dynamic evaporator model are Gaussian type denoted by $\mu$ in Figure 10 and Figure 11.

In addition, among different methods of defuzzifications, the centroid defuzzification method [42] is used in this model to convert the aggregated fuzzy set to a crisp output value $Y$ from the fuzzy set ${\Phi }^{\prime }$. This work computes the weighted average of the membership function or the center of gravity (COG) of the area bounded by the membership function curves [43]:

$$Y=\frac{{\int }_{}y.{\mu }_{{\Phi }^{\prime }}(y)dy}{{\int }_{}{\mu }_{{\Phi }^{\prime }}(y)dy}$$

(18)

The crisp value of $T_{r,o}$ and $T_{h,o}$ were calculated utilizing the above equation. The crisp values of the fuzzy inference systems are fed to the thermal inertia block to obtain the final output values of the fuzzy-based model.

The thermal inertia of the fluids in the evaporator is dependent on the fluid’s properties such as specific heat capacity, thermal conductivity, the volume flow rate of the fluid, etc., and the geometrical dimension and thermal properties of the evaporator wall. The thermal inertia of the fluid streams can be expressed by a time delay with respect to a change of the fluids temperature. The knowledge of the transient behavior of the evaporator model in the dynamic condition is used to define the time delay for the refrigerant and heat source temperature in the fuzzy-based model. The following transfer functions represent the time delay of the model.

Transfer function for the refrigerant’s outlet temperature delay is as follows:

$$\frac{0.999}{31s+1.0007}$$

(19)

The transfer function for the heat source outlet temperature delay is as follows:

$$\frac{1}{62s+0.99}$$

(20)

The numerator and denominator values of the two transfer functions were adjusted according to the response of the evaporator against the transient heat inputs studied in this research.

3.1.3. Evaporator Model Validation

The fuzzy-based dynamic model is validated with the data available from the dynamic FV model of the evaporator. The outputs of the fuzzy model with respect to the FV model are shown in Figure 14 and Figure 15. It can be observed from these figures that the fuzzy-based model can predict the evaporator outputs effectively in the dynamic condition. The data fitness values of 90.32% and 91.24% are obtained for the outputs of $T_{\_r\_o}$ and $T_{\_h\_o}$, respectively. The RMSE (Root Mean Squared Error) and MAPE (Mean Absolute Percentage Error) of the fuzzy model outputs with respect to the validation data is shown in

Table 5. Performance of fuzzy-based dynamic evaporator model.

Model Output	Fitness (%)	RMSE	MAPE (%)	Simulation Time (s)
				Fuzzy Model	FV Model
$T_{\_r\_o}$	90.32	1.10 (K)	0.19	5.19 s	3820.6 s
$T_{\_h\_o}$	91.24	3.09 (K)	0.58

. The fuzzy model can predict the outputs with a RMSE of 1.10 K for $T_{\_r\_o}$ and 3.09 K for $T_{\_h\_o}$, which gives a MAPE of 0.19% and 0.58%, respectively. Moreover, the time used for the simulation of the random heat source in the fuzzy-based model is 5.19 s, compared to 3820.6 s of the FV model. This implies that the fuzzy model with an accuracy of above 90% is 736.15 times faster than that of the FV model. These validation numbers indicate that the fuzzy logic concept is suitable to develop the dynamic evaporator model in the WHR system.

3.2. Other Component Models

3.2.1. Pump Model

The pump used in this research is a volumetric diaphragm pump [44], which is a positive displacement machine. The response of such pumps are much faster than the response of the heat exchanger in the ORC-WHR system [45]. For this reason, complicated dynamics of the pump are not necessary for the simulation of the WHR dynamics. A constant efficiency and the performance curve of the pump provided by the manufacturer can be used for the development of the pump model. The characteristic of the diaphragm pump is that the mass flow rate is proportional to the speed of the pump [14,46]. The performance curve of the pump can be used to derive the relationship between the speed of the pump and the refrigerant mass flow rate as follows [38,44]:

$$\frac{{\dot{m}}_{p1}}{{\dot{m}}_{p2}}=\frac{N_{p1}}{N_{p2}}$$

(21)

where ${\dot{m}}_p$ is the mass flow rate of the pump in kg/s and $N_p$ is the corresponding pump speed in RPM.

The specific work input as a function of the pump’s pressure and enthalpy difference can be calculated as follows:

$$\frac{W_p}{{\dot{m}}_p}=H_{p,o}^{}-H_{p,i}^{}=\frac{{\overline{\upsilon }}_p(P_{p,o}^{}-P_{p,i}^{})}{{\eta }_p}$$

(22)

where ${\overline{\upsilon }}_p$ is the specific volume in m³/kg; $H_{p,i}^{}$ and $H_{p,o}^{}$ are the pump inlet and outlet enthalpy in kJ/kg, respectively; $W_p$ is the pump work in kW; ${\eta }_p$ is the mechanical efficiency of the high pressure diaphragm pump; $P_{p,i}^{}$ and $P_{p,o}^{}$ are the inlet and outlet pressures of the pump in kPa, respectively.

3.2.2. Expander Model

An expander in a WHR system is generally one of two types: volumetric expander i.e., piston, scroll, etc. and a turbo-machinery type i.e., turbines. The time response of both volumetric and turbo-machinery expanders are much faster compared with the evaporator in the WHR system. Therefore, similarly with the pump, the dynamics of the expander is unnecessary in the simulation of WHR system. The steady state thermodynamic model of the expander is used for the simulation of the expander in this investigation. The specific work output of the expander is calculated as follows:

$$\frac{W_{\exp }}{{\dot{m}}_{\exp }}={\eta }_{\exp }(H_{\exp ,i}-H_{\exp ,o})$$

(23)

where $W_{\exp }$ is the expander work output, ${\dot{m}}_{\exp }$ is the refrigerant mass flow rate through the expander; $H_{\exp ,i}$ and $H_{\exp ,o}$ are the expander inlet and outlet enthalpy, respectively; ${\eta }_{\exp }$ is the expander efficiency.

3.2.3. Condenser Model

The dynamics of the WHR system are mostly formed by the dynamic behaviors of both evaporator and condenser. Therefore, the dynamic evaporator and condenser models were taken into account in the simulation of the WHR dynamics as shown in several reports [30,38,45]. Although both the evaporator and condenser have significant effects on the operating parameters of the ORC-WHR system, the dynamic model of the evaporator is of greater importance than that of the condenser. This is, especially at supercritical pressure, because of the unpredictable heat transfer mechanism at the evaporator. Therefore, a thermodynamic model of the condenser based on the state enthalpy is adequate to represent the condenser as follows:

$$\frac{Q_{con}}{{\dot{m}}_{con}}=(H_{\exp ,o}-H_{con,o})$$

(24)

where $Q_{con}$ is the condenser cooling power in kW, ${\dot{m}}_{con}$ is the mass flow rate of refrigerant through the condenser and $H_{con,o}^{}$ is the enthalpy at the condenser outlet.

3.2.4. Accumulator Model

To absorb the organic fluid fluctuations and make sure that there is liquid at the inlet of the pump, an accumulator is used between the pump inlet and the condenser outlet. The organic fluid level in the cycle is represented by following conservation equations [38,47].

$$\frac{d\phi }{dt}-\frac{{\dot{m}}_{ac,i}-{\dot{m}}_{ac,o}}{{\rho }_{r,l}{V}_{ac}}=0$$

(25)

$$\frac{d{H}_{ac}}{dt}-\frac{{\dot{m}}_{ac,o}(H_{ac,i}-H_{ac,o})}{{\rho }_{r,l}{V}_{ac}}=0$$

(26)

where $\phi$ is the relative liquid refrigerant level in the accumulator, ${\rho }_{r,l}$ is the liquid refrigerant density, $V_{ac}$ is the volume of the accumulator tank, ${\dot{m}}_{ac,i}$ is the mass flow rate of the refrigerant at the accumulator inlet, and ${\dot{m}}_{ac,o}$ is the mass flow rate of the refrigerant at the accumulator outlet. $H_{ac,i}$ and $H_{ac,o}$ are the enthalpy of the refrigerant at the inlet and outlet of the accumulator, respectively.

3.3. Model Integration

The component models of the supercritical ORC-WHR system described in Section 3.1 and Section 3.2 are integrated to form an overall system dynamic model. The overall model parameters are given in

Table 6. The organic Rankine cycle (ORC) model parameters.

Parameter	Symbol	Value	References
Pump efficiency	${\eta }_p$	0.75	[2,21]
Expander efficiency	${\eta }_{\exp }$	0.8	[2,38]
Evaporator heat transfer area	A	5.78 m2	[28]
Evaporator plate length	L	0.478 m
Evaporator plate width	W	0.124 m
Number of plates	Npl	100
Thermal conductivity	K	15 W/m K

. The inputs and outputs of each component in the dynamic model are linked according to the rationalities described as follows:

• The pump delivers the refrigerant mass flow rate at a proportion to the speed of the pump. Given the inlet and outlet pressure, the electrical power requirement to drive the pump and the temperature of the fluid at the outlet can be calculated.
• There is no enthalpy loss in between the pump and the evaporator. The heat recovery in the evaporator is a function of the inlet flow conditions and working pressure of the fluids.
• The expander can rotate freely without imposing any speed constraint. The amount of work output is a function of the enthalpy at the inlet and outlet of the expander.
• Provided the condenser outlet temperature is constant, the model calculates the required cooling power to achieve the desired temperature at the outlet.
• Given the inlet conditions and mass of the fluid, the accumulator maintains the outlet enthalpy and fluid level.

The pressure loss because of piping work between the evaporator and the expander in the overall WHR system is modelled using the Darcy–Weisbach pressure drop correlation as follows:

$$\Delta P_{pipe}=\frac{f_D\rho L_p{v}^2}{2D}$$

(27)

where $f_D$ is the Darcy friction factor, $\rho$ is the density of the pipe, $L_p$ is the length of the pipe, $v$ is the velocity of the fluid, and $D$ is the hydraulic diameter of the pipe.

The friction factor $f_D$ is calculated using the Haaland equation [48] as follows:

$$\frac{1}{\sqrt{f_D}}=-1.8\log \left[{\left(\frac{{\in }_p/D}{3.7}\right)}^{1.11}+\frac{6.9}{\mathrm{Re}}\right]$$

(28)

where ${\in }_p$ is the absolute roughness, and ${\in }_p/D$ is the relative roughness of the pipe.

4. Results and Analysis

The performance of the fuzzy-based overall model compared to the FV-based overall model in predicting the pressure at the expander inlet, the expander work output and the cycle efficiency is presented in this section. The simulation was conducted with the same input conditions as used for the fuzzy model validation in Section 3.1.3.

The pressure at the expander inlet after considering the loss is shown in Figure 16. It can be seen from this figure that the expander inlet pressure for both cases of the overall model are similar except for some small deviations around time 200 s, 400 s, and 800 s. These deviations of pressure are caused by the predicting error of the evaporator outlet temperature in the fuzzy domain, which is one of the functions of the pressure loss correlation used in Equation (27). However, the error in the fuzzy-based evaporator model predicting the expander inlet pressure was only 0.001%, which is shown in

Table 7. Mean values of operating parameters with respect to random inputs.

Parameters	FV Based Overall Model	Fuzzy Based Overall Model	MAPE (%)
$P_{\exp ,in}$	5996.64 kPa	5996.7 kPa	0.001%
$W_{\exp }$	2.53 kW	2.61 kW	3.11%
${\eta }_{cy}$	9.78	9.42	3.82%

. The performance of the overall models for the prediction of expander work output is shown in Figure 17. The MAPE of the fuzzy-based overall model for the work output is at 3.11% (see Table 7). This discrepancy between the fuzzy-based and FV-based overall model is due to the fact the output from the evaporator is also one of the inputs of the expander, thus the error associated in the evaporator model will lead to the error in the prediction of expander output.

The cycle efficiency of the ORC-WHR system is represented by ${\eta }_{cy}$ and can be written as:

$${\eta }_{cy}={\displaystyle \underset{t_1}{\overset{t_n}{\int }}\frac{W_{Net}}{Q_{ev}}}dt$$

(29)

where $W_{Net}$ is the net output work and $Q_{ev}$ is heat recovered at the evaporator.

The evaporator heat input $Q_{ev}$ in Equation (29) can be calculated from the mass flow rate and temperature of the refrigerant at the inlet and outlet of the fuzzy-based evaporator model.

The cycle efficiency of the fuzzy-based overall model in comparison to the FV-based model for the entire operation period is shown in Figure 18. Since the simulation in this paper was carried out with a random heat source and refrigerant input to the system, the efficiency was varied according to the combination of the heat source inputs, refrigerant inputs, and expander output, as shown in Figure 18. The cycle efficiency of the ORC-WHR system presented in this paper was neither controlled nor optimized during the simulation. Therefore, the outputs are followed according to the inputs only. It can be noticed that the maximum efficiency obtained in the simulation of the FV-based and fuzzy-based overall models are 13.36% and 12.86%, respectively. The FV-based overall model led to an average cycle efficiency of 9.78%, while the fuzzy-based overall model led to 9.42%. Since the cycle efficiency is a function of the temperature at the evaporator outlet, the error in the temperature calculation of the fuzzy-based evaporator model contributed to the efficiency calculation which can be seen in Figure 18 and Table 7. The results presented in this figure provide a MAPE error of less than 4% in predicting the cycle efficiency of the system using the fuzzy-based overall model compared to the FV model. Moreover, it can also be seen from Figure 18 that the efficiency prediction difference between the FV-based and the fuzzy-based model is always less than around 0.5%. The results indicate that the fuzzy-based evaporator model can be integrated with other components to predict the efficiency of the ORC-WHR system in dynamic conditions.

5. Conclusions

This paper presents the modelling of a dynamic ORC-WHR system using a well-known FV technique and a novel fuzzy method. The performance analysis of the dynamic model regarding individual components and overall model is presented.

Since the evaporator was the focus of the research, the dynamic model of this component was initially developed using the conventional FV method. The detailed performance analysis of the model regarding the transient response was carried out. Results show that the model is sufficiently stable and responses to the transient inputs are satisfactory. This implies that the model was robustly built to be used in the dynamic scenario of the WHR system. However, because of the high computation time of this method, a second model of the evaporator based on the fuzzy inference system was developed. Useful information about the evaporator in the dynamic condition such as transient responses, input-output ranges, thermal inertia, etc. obtained from the FV model were used to develop the fuzzy model. The outputs of the fuzzy-based model were validated with that of the FV model. The results from the validation indicate that the fuzzy inference system can be used to predict the evaporator outputs in the dynamic situation.

A pump, an expander, a condenser, and an accumulator were added to the proposed dynamic model of the evaporator to provide a complete overall model of the dynamic ORC-WHR system. The performance analysis of the overall model with the integration of the fuzzy evaporator model shows that the model was able to predict the pressure at the inlet of the expander, the expander output, and the cycle efficiency accurately. This shows that the fuzzy-based evaporator model in the WHR system can be used to simulate the system performance and to develop control strategies under transient heat input conditions.