Modeling for Organic Rankine Cycle Waste Heat Recovery System Development

This chapter introduces the modeling of organic Rankine cycle waste heat recovery (ORCWHR) system. The main goal of this chapter is to give an overview of ORC-WHR system modeling, especially focus on the heat exchanger models due to its key role in the ORCWHR system development. Six heat exchanger models considered in this chapter includes static model, 0D dynamic model, 1D finite volume model, 1D moving boundary model, 2D & 3D model. Model complexity, accuracy, and computation time are analyzed for the six heat exchanger models. More importantly, the heat exchanger model selection is discussed based on different phase of ORC-WHR system development, which facilitates the development of ORC-WHR system, and reduces the system development cost. In addition, a full ORC-WHR system model is presented as a modeling example including heat exchanger model, expander model, valve model and pump model.


Introduction
As the virtual representation of the ORC-WHR system, ORC-WHR model is a great tool to reduce the cost and time of system development. The ORC-WHR model can be categorized based on the heat exchanger modeling method as shown in Figure 1. Based on whether the model includes Ordinary Differential Equations (ODEs) or not, the heat exchanger models can be classified as dynamic model (w/ODEs) and static model (w/o ODEs). Under the dynamic category, the models are further classified based on the model dimension (0D, 1D, 2D and 3D). Under the 1D category, the models can be classified based on the concept of modeling (moving boundary model and finite volume model).
Generally, static modeling method only considers the energy balance of heat exchanger (i.e. the heat release from the heat sources equals to the heat absorption by the working fluid) [1][2][3]. Due to the static models lacking ODEs, there is no time varying parameter. Thus, this method can only analyze the steady state condition.
On the contrary, dynamic modeling methods are capable of simulating the transient conditions and provide the parameter vector that changes along the time vector. 0D model [4,5] lumps the parameter in a single point as shown in Figure 2 and the parameter is the same in any location of heat exchanger. For example, the heat source temperature is the same between the location near the exhaust gas inlet and the location near the exhaust gas outlet of heat exchanger. 1D model considers the one dimension in the direction of flow path. Two typical 1D models are finite volume model [6][7][8] and moving boundary model [9][10][11] as shown in Figures 3 and 4, respectively. For instance, the heat source temperature is different at different location of the flow path. The 2D and 3D models consider not only one dimension in the flow path direction, but also the directions perpendicular to the flow path axis.
The model is critical in the ORC-WHR system development. Static modeling method is utilized to analyze the energy flow and cycle efficiency at the beginning phase of the ORC system development. With the help of static models, heat sources selection, working fluid screening, expander machine selection, and cycle efficiency calculation can be roughly conducted. Dynamic models are developed in later phase of ORC-WHR system development to help  The ORC-WHR system model includes four main components and other supporting components. The four main components include evaporator, expander, condenser, and working fluid pump as shown in Figure 5. The other supporting components include valves, pipes, reservoir, feed pump, etc. The four main components exist in all the ORC-WHR system and the number of supporting components depend on the specific applications.
There are several challenges in the ORC-WHR system modeling:    The first modeling challenge is the heat exchanger modeling method selection. Heat exchanger is the key component of the system and there are several available modeling methods. The modeling process is time-consuming. Thus, it is extremely important to choose the right modeling method before diving into the modeling details.
The second challenge is the fidelity of the model or assumptions to be made. The more assumptions to be made, the less fidelity the model will be. On the other hand, the less assumptions require more modeling time and effort. Thus, there is a trade-off between the model fidelity and modeling effort. However, the smart choose of assumptions might significantly reduce the modeling time and effort depending on the purposes of the model.
The third challenge is the system model integration [7] and robustness. After the component models are finished, the component models need to be integrated together to simulate the entire ORC-WHR system. The inputs and outputs of the connected models must be compatible to each other. In addition, the robustness of the integrated model could be a problem in highly transient operating conditions such as ORC system warmup, cool down, or heat source fast step change, etc. Moreover, the coupling dynamics of working fluid temperature, pressure and phase change increases the difficulty of system model robustness improvement.
Section 2 gives an overview of the ORC-WHR system model and discusses the heat exchanger model comparison and selection at different phases of ORC-WHR system development. Section 3 presents the details of a finite volume heat exchanger model, expander model, valve model and pump model as an ORC-WHR system model example.

Overview of the ORC-WHR modeling
This section gives an overview of the ORC-WHR system modeling methods in terms of model complexity, accuracy, and computation time. With the modeling overview, it is easier to find the right modeling method in certain ORC-WHR system development phase. Among the four main components, pump model and expander model are generally simple compared with the heat exchanger model and there is no many choices respective to the pump and expander modeling methods. Heat exchanger model is the most challenging one in many cases. Thus, this chapter mainly focus on the heat exchanger modeling methods comparison and evaluation.

Heat exchanger model complexity
Model complexity indirectly represents the modeling time and effort required to build the model. The less complexity, the less modeling time and effort are required. In previous section, heat exchanger modeling methods are classified in Figure 1. In terms of model complexity, static heat exchanger model has the lowest complexity compared with dynamic model. The reason is that no ODE exists in the governing equations. All the equations are algebraic relation with energy balance. Using static method, the heat exchanger energy balance only has one equation and is given as follows: _ m hs c p, hs T hs, in À T hs, out where _ m hs is heat source mass flow rate, c p, hs is the heat source heat capacity, T hs, in =T hs, out are heat source heat exchanger inlet/outlet temperature, h f , in and h f , out are the working fluid enthalpy at heat exchanger inlet and outlet. They can be calculated based on the working fluid thermodynamic table as follows: where T f , in ,T f , out are working fluid heat exchanger inlet/outlet temperature, p f , evap is working fluid evaporation pressure. Different from static heat exchanger model, dynamic heat exchanger model considers the ODEs in the governing equations. In addition, the wall dynamics are included in the governing equations. Take 0D dynamic model as an example. Assuming there is not pressure drop across the heat exchanger, pressure dynamics are fast dynamics and can be neglected in the energy balance equation. The heat exchanger energy balance then include three equations: Heat source energy balance: m hs c p, hs dT hs dt ¼ _ m hs c p, hs T hs, in À T hs, out ðÞ À A hs, w U hs, w T hs À T w ðÞ (4) Working fluid energy balance: Wall energy balance: where A hs, w ,U hs, w are the heat transfer area and heat transfer coefficient between heat source and wall. The wall is the medium separating the heat source and working fluid. A f , w ,U f , w are the heat transfer area and heat transfer coefficient between working fluid and wall.
The 0D dynamic model has three equations in energy balance, whereas static model only has one equation in energy balance. Besides more equations in 0D dynamic model, each equation has more terms than that from static model. In addition, the working fluid mass balance equation is another equation in the 0D dynamic model as presented in Eq. (7). Thus 0D dynamic model is more complex than the static model and the 0D dynamic model requires more time and effort in modeling.
1D dynamic models share the same four governing equations (Eqs. (4-7)) with 0D in each discretized cell. The main difference is that 0D model has only one cell and 1D models have more than one cell. 1D dynamic models include finite volume model and moving boundary model. Finite volume model includes m discretized cells, each cell has the same volume.
Moving boundary model includes three cells, each cell has different volume, which are determined by the phase of working fluid. There are three phases in the working fluid inside the heat exchanger including pure liquid, mixed (liquid & vapor) and pure vapor. Each phase occupies one cell in the moving boundary model. The governing Eqs. (4)(5)(6)(7) are applied in each cell of 1D finite volume model and moving boundary model. Therefore, the 1D heat exchanger models has more equations than the 0D heat exchanger model. In terms of model complexity, 1D heat exchanger models are more complex than the 0D model. Even though both 0D and 1D models share similar governing equations, 1D models need to consider the boundary conditions between the adjacent cells, whereas 0D model only need to consider the boundary conditions at the heat exchanger inlet and outlet. Between the two 1D models, considering finite volume model generally has more than three cells, 1D moving boundary model has less governing equations than the 1D finite volume method and the number of different equations is (m-3)*4. Even though finite volume model has more equations, different cells share equations, which means as long as Eqs. 2D and 3D heat exchanger models are more complex than the 1D and 0D heat exchanger models and are generally modeled in CFD softwares (e.g. ANSYS, FLUENT, etc.).
Overall the heat exchanger model complexity rank can be given as follows: static model <0D dynamic model <1D dynamic finite volume model <1D dynamic moving boundary model <2D & 3D dynamic model.

Heat exchanger model accuracy
Accuracy is the model characteristic that everyone wants to maximize, because it determines the value of the model in some sense. Unlike model complexity, model accuracy can be easily quantified. There is a reference to compare with the model prediction and the accuracy represents the error between the model prediction and the reference.
Static heat exchanger model is utilized in the concept phase in the ORC-WHR development. In the concept phase, no components are selected and no experiments are conducted. Thus, in general, there is no reference data to evaluate the accuracy of the static model. Static model is usually utilized to assist basic energy balance between the heat sources and working fluid. This calculation only needs to give an estimate result and does not require high accurate model. In the static modeling process, many parameters are generally not considered such as heat transfer area and heat transfer coefficient in the heat exchanger, and heat loss from the heat exchanger to the ambient. Therefore, the static heat exchanger model accuracy is not high.
0D dynamic heat exchanger model makes less assumptions than the static model, which improves the model accuracy. To be specific, 0D model considers the heat transfer area and heat transfer coefficients in the heat exchanger. In addition, 0D model is generally validated by experimental data, which also increases the model accuracy compared with static model. 1D finite volume model and moving boundary model have less assumptions than the 0D model in that they consider one more dimension than the 0D model. With 1D model, the parameters of working fluid, heat source and wall have different values at different location in the axis of flow path, whereas 0D model share the same value at different location. This one more dimension feature equips the 1D models with higher fidelity than the 0D model. Thus 1D models have higher accuracy than the 0D model. Between the two 1D models, finite volume model has finer discretization resolution than the moving boundary model (m vs. 3 in Figures 3 and 4). Therefore, as the number m goes larger, the 1D finite volume model could have higher accuracy than the 1D moving boundary model.
2D and 3D heat exchanger model have less assumption than the 1D or 0D heat exchanger models, which equips them with higher accuracy.
Overall, the accuracy rank of all the heat exchanger models are as follows: Static model < 0D dynamic model < 1D moving boundary model < 1D finite volume model < 2D & 3D heat exchanger model.

Heat exchanger model computation time
The computation time is very important if the model needs to run online or the model is implemented in computational costly algorithms offline like Dynamic Programming. Static model only has algebraic equations and is the fastest model among all the heat exchanger models. The computation time of 0D and 1D models can be evaluated by the number of ODEs of the corresponding model. As mentioned earlier this section, the 0D model has the least number of ODEs, followed by 1D moving boundary model and 1D finite volume model, respectively. Similarly, the 2D model has more ODEs than 1D model and less ODEs than 3D model. Therefore, the computation time of all the heat exchanger models can be ranked as follows: static model <0D model <1D moving boundary model <1D finite volume model <2D & 3D model.

Model selection at different ORC-WHR system development phases
Selecting the right heat exchanger model at different phase of ORC-WHR system development has three benefits: 1. Meeting the certain development phase goals;

Reducing the time and effort of modeling;
3. Reducing the cost of modeling.
The ORC-WHR system development procedures can be explained using the diagram shown in Figure 6. The system development starts from concept design, in which phase the heat sources selection, working fluid selection, expander selection, expander power output form (electrical or mechanical) and system configuration are roughly evaluated and determined. It is common that some of the selection may be not finalized, which needs further investigation in the latter phases of the development. In the concept phase, some general energy balance are calculated to evaluate the power output at different operating condition and different system configuration. There is no experimental data and the calculation is not necessarily to be very accurate. Thus, the static heat exchanger model fits this development phase. The model is not very accurate, but its accuracy is enough to generate a general estimation of energy balance between heat sources and working fluid, and power output value from the expander machine.
After the concept phase, the component selection and development phase follows, during which the components hardware are selected based on the available products on the market or designed and manufactured. This is the first generation of the components selection, which may change several other generations based on the individual and integrated system experiments. Experimental data are collected from individual components. During this phase, the component models help the design of some components details such as the heat transfer area of heat exchanger, the nozzle area of turbine expander, the expansion ratio of piston expander, the size of the valve, the displacement of pump and the inner diameter of connected pipes. During this phase, the static model and 0D heat exchanger model help design the details of heat exchanger. The 0D model can be identified by the experimental data from the heat exchanger. The identified 0D model gives hints on updating the current heat exchanger generation to the next level by sweep some of key design parameters such as heat transfer area, section area of flow path, and length of flow path. The expander machine prototype can generate an efficiency or power map as a function of variables like expander speed, inlet pressure, outlet pressure, and expansion ratio. This map helps build either map-based or physics-based expander model. Figure 6. ORC-WHR system development procedures.

Organic Rankine Cycle Technology for Heat Recovery
The third development phase is system integration phase, during which the ORC-WHR components hardware are connected in the test rig. During this phase, the model is not required. As the components are integrated into an entire system, the next step is to conduct the experiments. However, without control, the ORC-WHR system experiments are hard to conduct due to the coupling of working fluid temperature, evaporation pressure, and the transient heat source.
The fourth development phase is the control development. The most important control in ORC-WHR system are the working fluid temperature control and evaporation pressure control. It is possible to develop the temperature and pressure control without a model (i.e. traditional PID feedback control). However, many simulation and experiments showed that the traditional PID feedback control cannot control the ORC-WHR system very well. Feedforward plus feedback or advanced controls (e.g. model predictive control) are proposed by many researchers in the field of ORC-WHR system. Both the feedforward and advanced controls require system models, either simple models or complex models. In the feedforward control, the highly accurate model helps improve the control performance. However, due to the combination with PID feedback control, the feedforward control model does not have to be very accurate. Because the feedback control helps correct the error brought by the feedforward control. Lower accuracy requirement helps reduce the modeling effort. Thus, static model or 0D model are common in feedforward control design. In advanced controls, the model accuracy has higher requirements than the feedforward control. It is because generally there is no feedback control to correct the model error. Accuracy is one of the constraints for advanced control and computation time is the other because the advanced control needs to run online. The accuracy constraint eliminates the possibility of selecting static model or 0D model and the computation time constraint eliminates the possibility of 1D finite volume model and 2D or 3D model. Thus, most of advanced controls developed for ORC-WHR system utilized 1D moving boundary model as the heat exchanger control-oriented model even though the moving boundary model has relative high model complexity.
The fifth development phase is the power optimization phase. After the control development, the experiments can be conducted. However, there is still a gap to reach the system design goal, which is the selection of the optimal reference trajectories for working fluid temperature at the outlet of evaporator heat exchanger, evaporation pressure, condensation pressure, working fluid subcooling temperature at condenser outlet etc. The model helps identify the optimal reference trajectories corresponding to the maximum expander power or net power from the ORC-WHR system at varying heat sources operating conditions. There are three methods to achieve the power maximization goal: (1) develop a map or a correlation between the optimal reference and the inputs such as heat sources mass flow rate and temperature. The map or correlation is from the optimal results from the steady state analysis with the help of model. In this method, accuracy is the most important factor and model complexity is also important. Therefore, finite volume method fits this requirement very well. (2) Develop an optimal reference trajectory based on transient driving cycle. Compared with steady states, transient driving cycle optimization requires more computation time, especially for dynamic programming algorithm. Thus, computation time is the most important factor for the model. The model accuracy should not be too low. Therefore, 0D model meets this criteria. (3) Directly optimize the power in the advanced control developed in the control development phase. In this case, no extra effort is needed in the power optimization phase. However, due to the computation time limitation of the online advanced control, the 'optimal' power calculated by the advanced control is local optimal rather than global optimal, which is the drawback of this third method.
As long as the control development and power maximization phases are done, the experimental implementation phase does not require the model.

Modeling of ORC-WHR system
This section presents the details of a full ORC-WHR system model, aiming to provide an example of modeling of entire ORC-WHR system. The configuration of the example system is shown in Figure 7. The tail pipe (TP) exhaust gas from an internal combustion engine is considered as heat source. Electrified turbine expander is chosen as expander machine. One valve is installed upstream of turbine expander to protect turbine from liquid working fluid during the system warmup or highly transient engine conditions. Another valve is installed in the bypass of turbine to allow working fluid bypass the turbine and also controls the working fluid evaporation pressure.

Heat exchanger modeling
Two heat exchangers exist in the ORC system including evaporator and condenser. In this chapter, only TP evaporator model is presented and the condenser is modeled using the same method. Two assumptions made in the TP evaporator model: (i) axial heat conduction are not considered in all three medium (working fluid, wall and TP exhaust gas), (ii) the wall temperature gradient in the radial direction is neglected, and (iii) pressure drop across the heat exchanger is not considered. The working fluid mass balance can be expressed as: where A cross represents cross-sectional area, subscript f represents working fluid, r represents density, _ m is mass flow rate, z represents spatial position in the axial direction. Mass flow balance in the exhaust gas side is ignored due to the exhaust gas fast dynamics. The working fluid and exhaust gas energy balance are expressed in the follow form: where p represent pressure, h represent enthalpy, d represents the effective flow path diameter, U represents the heat transfer coefficient, and ΔT represents temperature difference between the wall and the fluid (working fluid or exhaust gas).
The wall energy balance is shown below: where c p represents heat capacity, L represents the length in axial direction, A f , w represents the heat transfer area between working fluid and wall, U f , w represents the heat transfer coefficient between working fluid and wall. m η represents the heat exchanger efficiency multiplier, which accounts for heat loss to the ambient, A e, w represents the heat transfer area between exhaust gas and wall, U e, w is the heat transfer coefficient between exhaust gas and wall, and subscript w represents wall. Eqs. (8) and (9) are simplified to ODE Eqs. (4) and (7). Eq. (10), Eq. (4) and Eq. (7) are solved as follows: where k is the time step indices, Δt is length of time step. Overall, there are four equations to be solved for each cell: wall energy balance Eq. (11), working fluid mass balance Eq. (12), working fluid energy balance Eq. (13), and exhaust gas energy balance Eq. (13).
Heat transfer coefficients are calculated separately in working fluid side and exhaust gas side.
In the exhaust gas side, at each time step, heat transfer coefficient is calculated once and all the m cells share the same value, which is only a function of time. Eq. (14) is the expression of friction factor for the concentric tubes [12]:  (16) where ξ is friction factor, d in and d out are inner and outer diameters of concentric tube, respectively. The thermal conductivity of the exhaust gas is shown as follows: where d is hydraulic diameter, v d is dynamic viscosity, Pr is Prandtl number. Nusselt number expression, Eq. (20), of a concentric tube with insulated outer pipe wall is selected based on the heat exchanger structure [13]. where l is length of the pipe in the heat exchanger. The heat transfer coefficient between exhaust gas and wall are calculated with Eq. (21) [14]. The experimental evaporator construction differs slightly from concentric tubes, so a heat transfer coefficient multiplier (m U )i s applied.
U e, w ¼ m U Nu e k e d e For the working fluid side heat transfer coefficients, they are not only time dependent, but also space dependent due to the phase change of working fluid along the heat exchanger. The single phase heat transfer coefficients are calculated in Eq. (22). The expression is chosen based on geometry structure of the heat exchanger [13].
U f , w, i ¼