CET 96 DOI: 10.3303/CET2296030 Paper Received: 21 December 2021; Revised: 8 July 2022; Accepted: 26 August 2022 Please cite this article as: Muegge N., Kronberg A., Glushenkov M., Inguva V., Kenig E.Y., 2022, A Thermal Model for Recuperative Heat Engines Operating with Dense Working Fluids, Chemical Engineering Transactions, 96, 175-180 DOI:10.3303/CET2296030 CHEMICAL ENGINEERING TRANSACTIONS VOL. 96, 2022 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: David Bogle, Flavio Manenti, Piero Salatino Copyright © 2022, AIDIC Servizi S.r.l. ISBN 978-88-95608-95-2; ISSN 2283-9216 A Thermal Model for Recuperative Heat Engines Operating with Dense Working Fluids Nils Müggea, Alexander Kronbergb, Maxim Glushenkovb, Venkatesh Inguvaa, Eugeny Y. Keniga,* aChair of Fluid Process Engineering, Paderborn University, Pohlweg 55, 33098 Paderborn, Germany bEncontech B.V., TNW/SPT, P.O. Box 217, 7500 AE Enschede, The Netherlands eugeny.kenig@upb.de Currently, there is no commercially available technology for low temperature (below 100°C) heat conversion into mechanical or electrical energy. The Organic Rankine Cycle (ORC) is a typical technology for low temperature applications. However, ORC plants are rather complex and not economically sound for small-scale applications (<500kW). The isobaric expansion (IE) engine is a promising alternative to ORC plants for low- grade and medium-grade (>40°C) heat conversion. The IE engine works under external heat supply, and thus, almost every type of heat source is accessible. In this work, a mathematical model for transient thermal simulation of recuperative displacer-type IE engines operating with dense working fluids is proposed. The model is based on Nusselt number correlations for the heat transfer prediction and is coupled to the REFPROP 9.1 database yielding thermophysical properties. The first results of a parametric study performed with two different working fluids, namely R134a and carbon dioxide, are presented to show the impact of different operating frequencies on the heat transfer characteristics. 1. Introduction The version of an isobaric expansion (IE) engine proposed by Glushenkov et al. (2018) offers a promising potential for low and medium grade (between 40°C and 300°C) heat conversion into mechanical energy or electricity. This IE engine uses external heat supply and can be driven by any type of heat source, such as geothermal energy or industrial waste heat. Low temperature (<100°C) applications are particularly promising considering that the alternative Organic Rankine Cycle technology does not operate economically within this temperature region, mostly because of high installation costs. A key feature of the IE engine is the use of so- called dense working fluids, i.e. fluids that are fully liquid at low process temperature and gaseous or supercritical at high process temperature. The phase change effectively contributes to the engine power output, while most moving parts of the engine are located in regions filled with liquid, which reduces problems of sealing and lubrication. Glushenkov et al. (2018) proposed a special configuration comprising two displacer-type (or Bush-type) IE engines, which are interconnected by a recuperator as shown in Figure 1. This combination is called duplex arrangement. According to Figure 1a, the displacer in engine A moves towards the hot space (Hs) and displaces hot working fluid through the heat exchanger sequence comprising heater (H), recuperator (R) and cooler (C). At the same time, in engine B fluid is heated and receives heat from engine A in the recuperator. The flow direction of the fluid is indicated by small arrows. After one half cycle, the situation changes, and engine A receives heat from engine B, as depicted in Figure 1b. The thermodynamic behaviour of each IE engine can be described by an ideal rectangular cycle in the pressure-volume diagram. When working fluid is displaced from Cs to Hs, it is heated and the pressure inside the engine rises. As the pressure reaches a certain level 𝑝ℎ𝑖𝑔ℎ, the power piston starts to move and the overall internal volume of the engine expands, while the pressure inside remains constant. After having reached the final right position (volume occupied by Hs is at its maximum), the displacer begins its movement in the opposite direction while working fluid is displaced from Hs to Cs. The 175 working fluid is cooled, and hence, the pressure decreases until the pressure level 𝑝𝑙𝑜𝑤 is reached. Finally, the power piston moves back to its initial position due to continuous cooling of the working fluid at constant pressure. a) b) Figure 1: Schematic of two coupled IE engines in a duplex arrangement: engine A delivers heat to engine B during the first half cycle (a) and reversal of the flow direction in the second half cycle (b). The importance of well-designed heat exchangers in IE engine application has been emphasised by Knoke et al. (2020). They showed that dead volumes encountered due to internal volume of the heat exchanger equipment have a strong negative influence on the thermal efficiency of displacer-type IE engines. Therefore, they proposed using compact heat exchangers with high surface densities, such as the so-called printed circuit heat exchangers (PCHE). PCHEs consist of thin plates with semi-circular shaped channels, which are chemically etched into the plates as illustrated in Figure 2. Multiple plates are diffusion bonded in an alternating order, e.g., one plate with fluid delivering heat followed by a plate with fluid receiving heat. With this technology, compact heat exchangers with small diameter channels up to 0.5mm are feasible. Additionally, due to the diffusion bonding, these types of heat exchangers are highly pressure resistant which makes them useful for IE engine applications. In heat engine modelling, heat exchangers are commonly described based on idealised assumptions. Most often, ideal heat transfer conditions or a certain degree of recuperation are assumed, like in Invernizzi (2010). Such assumptions can lead to overestimation of the heat engine performance. Knoke et al. (2020) have proposed a steady-state design model for IE engine heat exchanger equipment using average values for the mass flow rates. However, their method is not capable of describing the strong dynamic behaviour of the cyclically phase-changing working fluid. Thus, the description of the heat transfer process must be deepened to provide a more accurate view on the overall performance of heat engines and especially the IE engine. In this work, a new model describing the dynamic heat transfer in compact heat exchangers within an IE engine application is presented. The method is used to provide information on the heat transfer for given heat exchanger designs and operating conditions. With this model, studies of the operating frequency influence on the heat transfer for subcritical refrigerant R134a and supercritical carbon dioxide are performed. 2. Mathematical model In order to investigate the influence of operating conditions on the transient heat transfer, a one-dimensional heat transfer model was developed. The model considers PCHEs for each of the heat exchangers H, R, and C (see Figure 1) and is based on the following assumptions for the heat exchanger sequence: • Heat transfer takes place only during the isobaric stages. • A recuperative IE engine setup is studied. The recuperator acts as a counter-flow heat exchanger, which connects engine A and engine B operating in counter-phase. • The temperature of the heating/cooling agent remains constant along the heater and cooler length, respectively. • The displacer moves cyclically between Hs and Cs. Its movement is approximated with a sinusoidal function. • Heat transfer between two fluid streams can be described based on overall heat transfer coefficients. 2.1 Treatment of printed circuit heat exchangers The channels of PCHEs have a semi-circular shape with diameter 𝐷; they are chemically etched into plates of thickness 𝐵, as shown in the cross-sectional view in Figure 2a. In the present work, heat transfer in Engine B Engine A H R C Hs Cs displacer power piston cylinder Engine B Engine A H R C 176 𝑥-direction is assumed to be negligible because the fluid temperature in parallel channels is nearly the same and the spacing between the channels is small. Based on this assumption, original channels are modelled as rectangular ones, as depicted in Figure 2b. The width of the rectangular channels is equal to the diameter of the original semi-circular geometry. This ensures the same heat transfer area in 𝑦-direction. The height 𝐻 of rectangular geometry is chosen in such a way that the cross-sectional area of a rectangular channel is equal to the area of the original channel, with the same flow velocity value. The rectangular geometry allows using geometrical periodicity shown in Figure 2b by the dash line. A more detailed view of one periodic element is given in Figure 2c. a) b) c) Figure 2: Schematic cross-sectional view on the original PCHE with semi-circular channels (a), simplified rectangular channels with a periodic element marked by the dash line (b) and the periodic element (c). 2.2 Governing equations It is assumed that the main part of heat transfer takes place during the isobaric stages. Consequently, the IE cycle is subdivided into a heating phase and a cooling phase. The pressure inside the engine is at a high pressure level during the heating phase and at a low pressure level during the cooling phase. For engine A starting with a cooling phase, the heating and cooling phases can be described as follows ℎ𝑒𝑎𝑡𝑖𝑛𝑔 𝑝ℎ𝑎𝑠𝑒: (𝛾 − 1)𝜏 ≤ 𝑡 < (𝛾 − 0.5)𝜏 𝑐𝑜𝑜𝑙𝑖𝑛𝑔 𝑝ℎ𝑎𝑠𝑒: (𝛾 − 0.5)𝜏 ≤ 𝑡 < 𝛾 𝜏 (1) where 𝛾 is the current cycle number (starting with 𝛾=1), 𝑡 is time and 𝜏 is the cycle period. Engine B is in a cooling phase when engine A is in a heating phase and vice versa. The model is presented for one of the engines of the duplex arrangement, while the description of the second engine is similar. The computational domain shown in Figure 3 is based on the periodic element depicted in Figure 2c and extended in 𝑧-direction. The following set of equations describes the working fluid inside the heat exchanger sequence. The continuity equation for the two-phase flow is 𝜕𝜌𝑡𝑝 𝜕𝑡 + 𝜕𝜌𝑡𝑝𝑤𝑧 𝜕𝑧 = 0 (2) in which the density of the two-phase flow is determined as 𝜌𝑡𝑝 = 𝜑𝜌𝐺 + (1 − 𝜑)𝜌𝐿 (3) where 𝜑 is the volume fraction of the gas phase and 𝜌𝐺 and 𝜌𝐿 are densities of the gas and liquid phase, respectively. The volume fraction can be recalculated from the vapour quality 𝑞 using the phase density values. Assuming that thermal conduction in flow direction is negligible compared to the convective energy transport, the energy transport equation is given by 𝜕𝜌𝑡𝑝𝑢 𝜕𝑡 + 𝜕𝜌𝑡𝑝𝑤𝑧ℎ 𝜕𝑧 = 𝑠𝑢 (4) in which 𝑢 is the specific internal energy, 𝑤𝑧 is the velocity in 𝑧-direction, ℎ is the specific enthalpy and 𝑠𝑢 is the energy source term due to heat transfer between two neighbouring channels in 𝑦-direction. Figure 3: Computational domain for one representative element in the flow direction (𝑧-axis). The specific enthalpy ℎ, fluid temperature 𝑇, vapour quality 𝑞 and phase densities 𝜌𝐺 and 𝜌𝐿 are obtained from the REPROP 9.1 database. 177 At the beginning of the simulation, Hs is at its maximum volume and the cooling phase is about to start. Consequently, the fluid is in high pressure state and the total mass of the working fluid that is displaced per half- cycle is 𝑚𝑠𝑤 = 𝑉𝑠𝑤 𝜌𝐺 (𝑇𝐻 , 𝑝ℎ𝑖𝑔ℎ) (5) where 𝑉𝑠𝑤 is the corresponding swept volume. It is assumed that the swept volume is twice as high as the internal volume of the heat exchanger sequence. This assumption is realistic according to Knoke et al. (2020). The volume change of Hs caused by the displacer movement leads to a mass flow entering the heater from the Hs or leaving it. Thus, the mass flow rate boundary condition is determined by the displacer position: �̇�(𝑧𝐻𝑠 , 𝑡) = 𝑚𝑠𝑤 𝑁𝐻 𝜋𝑓 sin(2𝜋𝑓𝑡) (6) where 𝑓 is the frequency and 𝑁𝐻 is the number of channels of the heater. The specific enthalpy of the fluid entering or leaving the domain is ℎ(𝑧𝐻𝑠 ) = { ℎ(𝑇𝐻 , 𝑝) �̇�(𝑧𝐻𝑠) > 0 ℎ(𝑇(𝑧𝐻 ), 𝑝) �̇�(𝑧𝐻𝑠) ≤ 0 (7) at Hs and ℎ(𝑧𝐶𝑠 ) = { ℎ(𝑇𝐶 , 𝑝) �̇�(𝑧𝐶𝑠 ) > 0 ℎ(𝑇(𝑧𝐶 ), 𝑝) �̇�(𝑧𝐶𝑠) ≤ 0 (8) at Cs, where 𝑧𝐻𝑠 is the coordinate at the interface between the heater and the hot space, 𝑧𝐶𝑠 is the coordinate at the interface between the cooler and the cold space, and 𝑧𝐻 and 𝑧𝐶 are the coordinates of the control volumes adjacent to Hs and Cs. The boundary conditions are given in Figure 3. The volumetric heat source term in Eq(4) was modelled based on a local heat balance using overall heat transfer coefficients. The overall heat transfer coefficient was defined for a periodic element in Figure 2c, assuming that the heat exchanger is made of steel. For channels filled with heating and cooling agent, constant heat transfer coefficients 𝛼𝐻/𝐶 =500W/(m²K) were assumed. For channels filled with working fluid, the heat transfer coefficients were obtained by Nusselt number correlations. In the present work, the correlations were primarily chosen based on recommendations for micro- channels and mini-channels by Kim and Mudawar (2014) given in Table 1. Table 1: Nusselt number correlations used in this work. Single-phase flow Condensing flow Flow-boiling Supercritical flow Gnielinski (1995) Akers and Rosson (1960) Gungor and Winterton (1986) Mokry et al. (2011) At the beginning of a heating or a cooling phase, the working fluid in the heater and cooler were initialised at heat source and heat sink temperature, respectively. In the recuperator, the temperature was initially set to the mean temperature between the heat source and heat sink. The velocity was set to zero in all segments. The model was implemented in MATLAB®. 3. Results and Discussion Cyclic heating and cooling of subcritical R134a and supercritical carbon dioxide was investigated for four different operating frequencies varied between 0.125Hz and 1Hz. The geometric specifications of the heat exchanger equipment are listed in Table 2; they were kept for all simulations. Along with the channel diameter 𝐷 and the plate thickness 𝐵, each heat exchanger of H, R, C is characterised by its length 𝐿 and total number of channels per fluid stream 𝑁. The analysis was carried out for a heat source temperature of 𝑇𝐻=353K and a heat sink temperature of 𝑇𝐶 =283K. Table 2: Geometric specification of the heat exchangers. 𝐷/mm 𝐵/mm 𝐿/m 𝑁 Heater (H) 0.5 2 0.3 400 Recuperator (R) 0.5 2 0.3 400 Cooler (C) 0.5 2 0.3 500 The low pressure during the cooling phase with working fluid R134a was 𝑝𝑙𝑜𝑤=5.7bar and the high pressure during the heating phase was 𝑝ℎ𝑖𝑔ℎ=24.2bar. The pressure levels were chosen based on the saturation temperature of the working fluid. At low pressure, the saturation temperature was 10K above 𝑇𝐶 and at high pressure, the saturation temperature was 10K below 𝑇𝐻. This ensures a temperature difference to drive the heat 178 transfer process during the phase change. In the study with carbon dioxide, the low pressure was 65bar and the high pressure was 100bar. A mesh independence study was carried out yielding the value 1.5mm/segment. The simulation results were used to evaluate the heat exchanger efficiency, which is defined as 𝜖𝐶𝑃 = 100 ⋅ ℎ(𝑇𝐻,𝑝𝑙𝑜𝑤 )−ℎ(𝑧𝐶) ℎ(𝑇𝐻,𝑝𝑙𝑜𝑤 )−ℎ(𝑇𝐶 ,𝑝𝑙𝑜𝑤 ) . (9) for the cooling phase (CP). The evaluation of the heating phase (HP) is similar. The heat exchanger efficiency for one half-cycle with working fluid R134a operating at 0.125Hz is depicted in Figure 4a. The solid blue line shows the efficiency of CP (engine A) while the efficiency during the HP (engine B) is given by the solid red line. The efficiency of the recuperation is represented by the blue and red dash lines for CP and HP, respectively. Initially, the efficiency of both HP and CP remains at 100% for about one second, because the heater and cooler are initialised at heating/cooling agent temperatures and displacer movement is slow. Afterwards, convection increases with higher displacer speed, bringing cold fluid to the hot side, which leads to a decreased HP efficiency in engine B. The CP efficiency of engine A reacts accordingly. At 1.3s, the HP is at its lowest efficiency, because of the reduced residence time. At the end of the half-cycle, both CP and HP efficiencies approach 100%. Figure 4b shows the heat exchanger efficiencies of the HP (red) and CP (blue) at four different operating frequencies. The efficiencies were evaluated at the end of a half-cycle. The highest efficiencies are reached at the lowest frequency of 1/8Hz, while the lowest efficiencies are observed for the fastest displacer movement at 1Hz. It is visible that the reduction of the HP efficiency, dropping from 98% to 58%, is lower than that of CP, decreasing from 99% to 36%. a) b) Figure 4: Heat exchanger efficiency for working fluid R134a operating at 0.125Hz (a) and the heat exchanger efficiencies at the end of a half-cycle for four different frequencies (b). a) b) 179 Figure 5: Heat exchanger efficiency for carbon dioxide operating at 0.125Hz (a) and the heat exchanger efficiencies for four different frequencies (b). Figure 5a depicts the heat exchanger efficiency of a recuperative half-cycle with working fluid carbon dioxide operating at 0.125Hz. In contrast to the refrigerant R134a, the high pressure level in this case is above the fluid critical pressure. Both the CP (blue) and HP (red) efficiency arrive at their minimum after the maximum displacer speed is reached at 2s. While the displacer slows down in the interval between 2s and 4s, the efficiencies of CP and HP increase until they reach their final level at 97% and 99%, respectively. The efficiencies of the recuperation increase during the first 1.5s. This is attributed to the hot fluid entering the recuperator from the heater in engine A providing heat for the HP in engine B. Similar effects were obtained for R134a. At the end of the half-cycle, the recuperation efficiencies reach 36% for CP and 54% for HP, which are at a higher level than in the case of R134a. The heat exchanger efficiencies of CP and HP for four different frequencies are shown in Figure 5b. Generally, an efficiency decrease with rising frequencies is observed, although this decrease is not as significant as for R134a. In particular, the efficiency drop of the heating phase is lower, which is attributed to the enhanced heat transfer and better recuperation with carbon dioxide. 4. Conclusions A new thermal model for compact heat exchangers in IE engine application was developed. First investigations with two different working fluids, namely refrigerant R134a in subcritical state and supercritical carbon dioxide, were carried out to study the heat exchanger performance at four different frequencies between 0.125Hz and 1Hz. In general, high heat exchanger efficiencies are established for lower operating frequencies. This makes IE engines favourable for practical applications, because of their ability to have a high power density per cycle. The efficiency is more sensitive to increasing frequencies for working fluid R134a than for carbon dioxide. Furthermore, this efficiency is higher with carbon dioxide achieving values between 62% at 1Hz and 99% at 0.125Hz in the heating phase, which is favourable for optimal utilisation of the heat sources. To quantify the relationship between the frequency and the heat exchanger efficiency, further parametric studies are necessary. Nomenclature 𝐵 – thickness, m 𝐷 – diameter, m 𝑓 – frequency, Hz ℎ – specific enthalpy, J/kg 𝐻 – height, m 𝐿 – length, m �̇� – mass flow rate, kg/s 𝑚𝑠𝑤 – swept mass, kg 𝑁 – number of channels, - 𝑝𝑙𝑜𝑤 – low pressure level, Pa 𝑝ℎ𝑖𝑔ℎ – high pressure level, Pa 𝑞 – vapour quality, - 𝑠𝑢 – energy source term, W/m³ 𝑡 – time, s 𝑇 – working fluid temperature, K 𝑇𝐶 – heat sink temperature, K 𝑇𝐻 – heat source temperature, K 𝑢 – specific internal energy, J/kg 𝑉𝑠𝑤 – swept volume, m³ 𝑤𝑧 – velocity of the fluid phase, m/s 𝑥, 𝑦, 𝑧 – cartesian coordinates, m 𝛼 – heat transfer coefficient, W/(m2·K) 𝛾 – cycle number, - 𝜖 – heat exchanger efficiency, - 𝜌𝐺 – gas-phase density, kg/m³ 𝜌𝐿 – liquid-phase density, kg/m³ 𝜌𝑡𝑝 – two-phase flow density, kg/m³ 𝜏 – cycle period, s 𝜑 – vapour volume fraction, - Acknowledgement This project has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No 764042. References Akers W.W., Rosson H.F., 1960, Condensation inside a horizontal tube, Chemical Engineering Progress Symposium, 56, 145-149. Glushenkov M., Kronberg A., Knoke T., Kenig E.Y., 2018, Isobaric expansion engines: new opportunities in energy conversion for heat engines, pumps and compressors, Energies, 11, 154-175. Gnielinski V., 1995, A new mathematical calculation method for the heat transfer in the transition regime between laminar and turbulent pipe flow (in German), Forschung im Ingenieurwesen, 61, 240-248. Gungor K.E., Winterton R.H.S., 1986, A general correlation for flow boiling in tubes and annuli, International Journal of Heat and Mass Transfer, 29, 351-358. Invernizzi C.M., 2010, Stirling engines using working fluids with strong real gas effects, Applied Thermal Engineering, 30, 1703-1710. Kim S.-M., Mudawar I., 2014, Review of database and predictive methods for heat transfer in condensing and boiling mini/micro-channel flows, International Journal of Heat and Mass Transfer, 77, 627-652. Knoke T., Kronberg A., Glushenkov M., Kenig E.Y., 2020, On the design of heat exchanger equipment for novel- type isobaric expansion engines, Applied Thermal Engineering, 167, 114382. Mokry S., Pioro I., Farah A., King K., Gupta K., Peiman W., Kirillov P., 2011, Development of supercritical water heat-transfer correlation for vertical bare tubes, Nuclear Engineering and Design, 241, 1126-1136. NIST, Reference Fluid Thermodynamic and Transport Properties Database (REFPROP); Version 9.1, https://www.nist.gov/srd/reprop, 2019. 180 76muegge.pdf A Thermal Model for Recuperative Heat Engines Operating with Dense Working Fluids