Improvement of a Nusselt-Based Simulation Model for Heat Transfer in Rotary Heat Exchangers

: In the last 50 years, the technology of rotary heat exchangers has not changed considerably. A reliable simulation can help improve the design of this technology. In this work, a simulation for rotary heat exchangers was developed and validated with multiple experimental data. This simulation takes an innovative approach based on locally calculated heat transfer coe ﬃ cients and considers the entry region e ﬀ ect. This approach proved to be accurate since the average di ﬀ erence between the experimental results and the proposed model with a constant heat boundary condition is 0.1% and the maximum absolute deviation 1%. Experimental, as well as simulation results, indicate that lower empty tube gas velocity (1 m / s) and higher rotational speed (12 rpm) improve thermal e ﬃ ciency compared to commonly used operating conditions. Additionally, a new model for predicting the local internal Nusselt number for sine ducts in the rotor channels is proposed, which considers the entry region e ﬀ ect.


Introduction
Rotary heat exchangers are already among the most efficient [1,2] and compact [3] technologies for air-to-air heat recovery. Over the last decade, heat transfer efficiency requirements in the European Union (EU) for air-conditioning systems have been steadily increasing. In the EU after 2020 it will be legally required in non-residential buildings that the thermal efficiency must be greater than 85% [4]. The thermal efficiency η refers to the relationship between the actual heat flow transferred by the heat exchanger . Q trans f erred and the maximal heat flow that is theoretically possible in an adiabatic and infinite long heat exchanger . Q max as defined by different authors [5,6] η ≡ . Q trans f erred . Q max (1) This definition of efficiency does not include electricity and mechanical power losses due to ventilators in the system. These loses are out of the scope of this work and a possible point for research in further work.
The relevance of these legislative efficiency measures relies on the amount of installed devices. According to a study by Kaup and Kampeis [7], in Germany alone more than 60,000 ventilation systems were delivered in 2009 and approximately 60% of them had a heat recovery unit. By 2012, deliveries had doubled and more than 80% of the systems were equipped with a heat recovery unit. This development is due to higher energy costs, since ventilation represents between 35% and 38% of energy loss in buildings [7]. These numbers do not consider the demand for industrial purposes and yet show how significant a thermal efficiency improvement on heat recovery systems can influence energy demand.
Among the different types of heat recovery system, rotary heat exchangers can provide thermal efficiency greater than 80%, making them the most efficient technology known to date [8,9]. This technology was originally developed and patented in 1922 by the Swedish engineer Frederik Ljungström [10]. In the following decades, thermodynamic, aerodynamic and structural engineering principles were examined to further develop rotary heat exchangers [11]. Rotary heat exchanger innovations stagnated during the 1960s. Since then, for example the rotor length in air ventilation systems has been a fixed 20 cm as a standard design parameter. If a change in these design or operational parameters could increase the thermal efficiency by 5-10%, this would already meet the European target [5], making it possible for rotary heat exchangers to comply with the regulations in the future. This would open the possibility of saving more energy, and decrease costs and CO 2 emissions with regard to the current rotary heat exchangers.
The purpose of this work is to develop and validate a simulation of rotary heat exchangers that can predict their thermal efficiency by a given set of operational and design parameters. With this simulation, a better design of rotary heat exchangers can be achieved in the future. Previous work has tried to fulfil this task; the accuracy of these state-of-the-art simulations was in the best cases limited (maximum absolute deviation of 3.5% points for Leong [12] or 3.71% in average for Özdemir [13]). To achieve the goal of this work, a high-precision simulation of the heat transfer for rotary heat exchangers, with the following innovative strategies are implemented in this work: (1) To improve the accuracy of the simulation results of rotary heat exchangers by using a locally calculated Nusselt number model in order to determine the heat transfer coefficients between the gas streams and the solid material. To date, many studies have used constant heat transfer coefficients along the tube, where the heat transfer takes place [13][14][15]. (2) Since there is no literature on local Nusselt numbers for the developing region in sine ducts, this effect is usually neglected in other works. In this work, this knowledge gap is addressed. (3) For validating the developed simulation model, multiple experimental data will be obtained from a pilot plant and from a publicly available online database from Eurovent [16]. This will ensure that the developed simulation model is validated with data from different rotary heat exchangers. Previous simulations were only validated with one rotary heat exchanger [2,12,15,[17][18][19][20][21][22] or a single rotor from a literature source [13,[23][24][25]. (4) Previous authors have used constant thermophysical properties [12,13,15,19,20,23,26].
Here, for improving simulation accuracy, thermophysical properties are calculated locally and are temperature dependent. (5) The state-of-the-art simulations for a rotary heat exchanger use either the ε-NTU method [2,12,14,15,17,[20][21][22][23][24][27][28][29][30][31][32] or computer fluid dynamics [13,18,25,26]. Finite difference methods, with locally calculated Nusselt number models and entry-region considerations, as done in this work, have until now been non-existent in literature. In the current context, the development of a more accurate simulation model for rotary heat exchangers is highly relevant because on the one hand, it predicts the extent of the impact of each parameter on the operation parameters, in particular the heat recovery efficiency of the system. On the other hand, such a model can serve as a basis for optimization. Without this accuracy, a further optimization would be futile. This is of highest significance because it can satisfy both sustainable manufacturing of this system and compliance with regulations regarding energy recovery. These requirements have become more and more strict over the last decades and probably will do so in future too.

Materials and Methods
For practical reasons, rotary heat exchangers, which are also known in the literature as "thermal wheels", "heat recovery wheels", "Kyoto wheels" or simply "rotors" will in most cases hereafter be referred to as "rotor" or "thermal wheel".
Thermal wheels are placed inside housing that is divided into two sectors, as shown in Figure 1. The thermal wheel rotates continuously from one sector to the other. In one sector, hot process air (1.1) flows into the wheel, the wheel metal material is in contact with the hot process air, which transfers heat to the wheel, the wheel heats up and stores the energy provided by the process air. The process air flows continuously through the wheel leaving the wheel at a lower temperature (exhaust air 1.2). Since the wheel is continuously rotating, the heated part of the wheel rotates into the other sector where it meets the colder ambient air (2.1). The heated solid material transfers the stored heat to the ambient air. The ambient air leaves the wheel at a higher temperature (supply air 2.2). The now cold wheel section continues rotating to the first sector where it again heats up. Since the thermal wheel is rotating continuously from one sector into the other, one half of the wheel is continuously being heated and storing energy from the process air, while the other half is continuously cooling and returning the stored energy to the ambient air.

Materials and Methods
For practical reasons, rotary heat exchangers, which are also known in the literature as "thermal wheels", "heat recovery wheels", "Kyoto wheels" or simply "rotors" will in most cases hereafter be referred to as "rotor" or "thermal wheel".
Thermal wheels are placed inside housing that is divided into two sectors, as shown in Figure 1. The thermal wheel rotates continuously from one sector to the other. In one sector, hot process air (1.1) flows into the wheel, the wheel metal material is in contact with the hot process air, which transfers heat to the wheel, the wheel heats up and stores the energy provided by the process air. The process air flows continuously through the wheel leaving the wheel at a lower temperature (exhaust air 1.2). Since the wheel is continuously rotating, the heated part of the wheel rotates into the other sector where it meets the colder ambient air (2.1). The heated solid material transfers the stored heat to the ambient air. The ambient air leaves the wheel at a higher temperature (supply air 2.2). The now cold wheel section continues rotating to the first sector where it again heats up. Since the thermal wheel is rotating continuously from one sector into the other, one half of the wheel is continuously being heated and storing energy from the process air, while the other half is continuously cooling and returning the stored energy to the ambient air. is the energy that the rotor provides to the supply air, which is based on the enthalpy change from the ambient air (2.1) to the supply air flow (2.2). Its value is identical to the energy decrease of the process air (1.1) to the exhaust air (1.2). If the heat exchanger is considered without heat losses, the energy balance for the fluid streams can be expressed as: is defined as the energy that would be transferred if the heat exchanger were adiabatic and infinitely long.
is limited by the stream with the lower mass flow: Since both mass flows in this paper are assumed to be equal on each side and in narrow temperature ranges, the specific heat capacity can be considered constant. The equation is commonly presented as: The right side of Equation (4) is what is often called "temperature effectiveness" by Shah and Sekulić [6] or "temperature efficiency" by Eurovent [16] and for comparison reasons this will be the most important parameter to be further used in this work. Hereafter it is only referred to as "temperature effectiveness" . . Q trans f erred is the energy that the rotor provides to the supply air, which is based on the enthalpy change from the ambient air (2.1) to the supply air flow (2.2). Its value is identical to the energy decrease of the process air (1.1) to the exhaust air (1.2). If the heat exchanger is considered without heat losses, the energy balance for the fluid streams can be expressed as: .
. Q max is defined as the energy that would be transferred if the heat exchanger were adiabatic and infinitely long.
. Q max is limited by the stream with the lower mass flow: Since both mass flows in this paper are assumed to be equal on each side and in narrow temperature ranges, the specific heat capacity can be considered constant. The equation is commonly presented as: The right side of Equation (4) is what is often called "temperature effectiveness" by Shah and Sekulić [6] or "temperature efficiency" by Eurovent [16] and for comparison reasons this will be the most important parameter to be further used in this work. Hereafter it is only referred to as "temperature effectiveness" η.

Research Tools
The methodological approach in this paper uses simulations based on thermodynamic equations as stated in Section 2.1.1 and a pilot plant set-up as depicted in Section 2.1.2. These are be the main tools used for obtaining and analyzing data afterwards.

Simulation
In this section, the simulation tool and methods are explained: this comprises the energy balance (Section 2.1.1.1), assumptions and boundary conditions (Section 2.1.1.2), the simulation concept and numeric, as well as the thermodynamic and fluid dynamic considerations.

Energy Balance
The base model for the heat transfer of rotors has already been explained in detail in existing literature [30,33]. It is based on the energy balances between a metal channel and the gas that flows through it. This model assumes: • One-dimensional heat conductivity in axial direction z for the gas and solid phase (see Figure 2 for the z direction).

•
The solid temperature is homogeneous in the radial direction r as depicted on Figure 2 (infinite radial thermal conductivity). Please note that the radial direction is perpendicular to the z coordinate. • At each given cross-section inside the duct along z the gas is perfectly mixed and has a homogeneous temperature distribution.

•
The solid phase is homogeneous in thermophysical properties at the start.

•
No back mixing in the gas flow takes place.

•
The system is considered to have no heat losses. • Potential and kinetic energy effects are negligible.

•
Only one-dimensional gas flow in the direction z is considered.

•
The gas is incompressible and at constant pressure.

•
Enthalpy is only a function of temperature.

Research Tools
The methodological approach in this paper uses simulations based on thermodynamic equations as stated in Section 2.1.1 and a pilot plant set-up as depicted in Section 2.1.2. These are be the main tools used for obtaining and analyzing data afterwards.

Simulation
In this section, the simulation tool and methods are explained: this comprises the energy balance (Section 2.1.1.1), assumptions and boundary conditions (Section 2.1.1.2), the simulation concept and numeric, as well as the thermodynamic and fluid dynamic considerations.

Energy Balance
The base model for the heat transfer of rotors has already been explained in detail in existing literature [30,33]. It is based on the energy balances between a metal channel and the gas that flows through it. This model assumes: • One-dimensional heat conductivity in axial direction for the gas and solid phase (see Figure 2 for the direction).

•
The solid temperature is homogeneous in the radial direction as depicted on Figure  2 (infinite radial thermal conductivity). Please note that the radial direction is perpendicular to the coordinate. • At each given cross-section inside the duct along the gas is perfectly mixed and has a homogeneous temperature distribution.

•
The solid phase is homogeneous in thermophysical properties at the start.

•
No back mixing in the gas flow takes place.

•
The system is considered to have no heat losses.

•
Potential and kinetic energy effects are negligible.

•
Only one-dimensional gas flow in the direction is considered.

•
The gas is incompressible and at constant pressure.

•
Enthalpy is only a function of temperature.

•
No mass transfer takes place between gas and solid, (Please note that if . = . , then thermal efficiency and temperature effectiveness are equal).
• No adsorption, absorption or condensation is considered. The energy balances result in the following partial parabolic differential equation system which can be found easily in literature [34]: The energy balances result in the following partial parabolic differential equation system which can be found easily in literature [34]: Energies 2021, 14, 10 5 of 26

. Boundary Conditions
For the solid phase, no heat transfer takes place in the axial direction z at the entrance and at the outlet of the channel, following the same assumptions of Gnielinski [33] and Vohl [35]: Based on previous work, at the outlet for the gas phase, adiabatic conditions are assumed [36]: Based on previous literature [36], from the energy balance at the inlet for the gas phase, we obtain:

Design and Operational Parameters in the Simulation
Rotors are typically made of aluminum. Two bands of metal foils, typically of 20 cm width, one corrugated and one flat are coiled as shown in Figure 3. This honeycomb-similar structure has many small flow channels that allow gas to pass through in the axial direction.

Boundary Conditions
For the solid phase, no heat transfer takes place in the axial direction at the entrance and at the outlet of the channel, following the same assumptions of Gnielinski [33] and Vohl [35]: Based on previous work, at the outlet for the gas phase, adiabatic conditions are assumed [36]: Based on previous literature [36], from the energy balance at the inlet for the gas phase, we obtain:

Design and Operational Parameters in the Simulation
Rotors are typically made of aluminum. Two bands of metal foils, typically of 20 cm width, one corrugated and one flat are coiled as shown in Figure 3. This honeycomb-similar structure has many small flow channels that allow gas to pass through in the axial direction. As shown in Figure 3, due to the manufacturing process, the metal foils are coiled into a spiral. In a spiral the curvature along the radius decreases continuously. For example, the outer edge of the rotor has the lowest curvature. Due to this curvature effect, the wave crests (highest point of a wave) stretch apart and the wave troughs (lowest point of a wave) are compressed. Therefore, the distance between the wave troughs is smaller than the distance between the wave crests. Due to this effect, the wavelength is not equal at wave crest and wave trough. This is a challenge that in practice would make it unfeasible to develop a model considering each individual tube inside the rotor. Therefore, this work As shown in Figure 3, due to the manufacturing process, the metal foils are coiled into a spiral. In a spiral the curvature along the radius decreases continuously. For example, the outer edge of the rotor has the lowest curvature. Due to this curvature effect, the wave crests (highest point of a wave) stretch apart and the wave troughs (lowest point of a wave) are compressed. Therefore, the distance between the wave troughs is smaller than the distance between the wave crests. Due to this effect, the wavelength is not equal at wave crest and wave trough. This is a challenge that in practice would make it unfeasible to develop a model considering each individual tube inside the rotor. Therefore, this work simplifies the geometry as shown in Figure 4, which can be described based on the following design parameters: Material thickness s: thickness of the foil used to manufacture the rotor.

•
Wave height h * : distance between two coils in which the waves are formed. • Wave radius r: radius of the circle that is formed by the corrugated foil when the waved layer changes direction. • Inclination angle β: refers to the angle formed between the straight part of the wave and the flat aluminum layer. • Rotor length l: distance that the gas flows inside the rotor as shown in Figure 2. • Wavelength t * : the distance between two wave crests; dependent on the previous design parameters β, r, h * , s. Since the model assumes a flat geometry, distances at wave crest, wave through and wave middle are identical.

•
Wave height-to-length ratio ς is used for certain heat transfer models and is defined as follows: Energies 2020, 13, x FOR PEER REVIEW 6 of 26 simplifies the geometry as shown in Figure 4, which can be described based on the following design parameters: • Material thickness : thickness of the foil used to manufacture the rotor.

•
Wave height ℎ * : distance between two coils in which the waves are formed.

•
Wave radius : radius of the circle that is formed by the corrugated foil when the waved layer changes direction. • Inclination angle : refers to the angle formed between the straight part of the wave and the flat aluminum layer.

•
Rotor length : distance that the gas flows inside the rotor as shown in Figure 2. • Wavelength * : the distance between two wave crests; dependent on the previous design parameters , , ℎ * , . Since the model assumes a flat geometry, distances at wave crest, wave through and wave middle are identical.

•
Wave height-to-length ratio is used for certain heat transfer models and is defined as follows: Using the parameters described, the length of a waved duct * , which is within the length * , is calculated as follows: Based on this geometrical description, the void fraction and the volumetric specific surface is defined as follows: Operating parameters also have an impact on rotor thermal efficiency : • Rotational speed of the rotor . • Empty tube gas velocity is the gas velocity before the rotor inlet. Since the rotor partially covers the flow area, it is reducing the cross-section area as shown in Figure  2, the internal gas velocity is given by: Using the parameters described, the length of a waved duct l * , which is within the length t * , is calculated as follows: Based on this geometrical description, the void fraction ε and the volumetric specific surface a v is defined as follows: Operating parameters also have an impact on rotor thermal efficiency η: • Rotational speed of the rotor n. • Empty tube gas velocity v a is the gas velocity before the rotor inlet. Since the rotor partially covers the flow area, it is reducing the cross-section area as shown in Figure 2, the internal gas velocity v in is given by: The value of gas velocity v is affected by density changes caused by variations of temperature along the channel. Nevertheless, the mass flow remains constant. Therefore, in the present work, sometimes the velocity will be shown in mass flux (mass flowrate per surface area expressed in kg·m −2 ·s −1 ) as follows:

Simulation Concept and Numerical Methods
For the simulation, the following assumptions were made: • Equal and steady state gas mass flow on both sides of the rotor.

•
No leakages at the rotor sealing.

•
No purge sector considered.
The heat transfer mechanism of the whole rotor can be described based on the heat transfer in a single channel: Hot air is blown into the channel and heats up the channel for a certain time. This time span corresponds to the time of half a revolution of the rotor. The channel then enters the other sector of the housing where cold air is blown through the channel in the opposite direction for the same period. This cold air heats up by drawing the energy from the channel.
The process of cooling down and heating up of a channel is repeated continuously. The gas outlet temperature of a single channel is continuously changing along the rotation as shown in Figure 5. However, after the rotor system has reached steady-state, at any given angle the outlet temperature from a channel is constant over time. Therefore, the resulting outlet gas temperature of the stream composed by the individual flow of each channel remains constant, which is the average outlet gas temperature over a half-cycle in Figure 5.
The value of gas velocity is affected by density changes caused by variations of temperature along the channel. Nevertheless, the mass flow remains constant. Therefore, in the present work, sometimes the velocity will be shown in mass flux (mass flowrate per surface area expressed in kg·m −2 ·s −1 ) as follows:

. Simulation Concept and Numerical Methods
For the simulation, the following assumptions were made: • Equal and steady state gas mass flow on both sides of the rotor.

•
No leakages at the rotor sealing.

•
No purge sector considered.
The heat transfer mechanism of the whole rotor can be described based on the heat transfer in a single channel: Hot air is blown into the channel and heats up the channel for a certain time. This time span corresponds to the time of half a revolution of the rotor. The channel then enters the other sector of the housing where cold air is blown through the channel in the opposite direction for the same period. This cold air heats up by drawing the energy from the channel.
The process of cooling down and heating up of a channel is repeated continuously. The gas outlet temperature of a single channel is continuously changing along the rotation as shown in Figure 5. However, after the rotor system has reached steady-state, at any given angle the outlet temperature from a channel is constant over time. Therefore, the resulting outlet gas temperature of the stream composed by the individual flow of each channel remains constant, which is the average outlet gas temperature over a half-cycle in Figure 5.  Figure 5 shows that where the rotor enters the cold sector (for angles in the range 0° to 5°), the temperature outlet profile shows a higher gradient. The reason for this is the carry-over gas from the hot sector. The time it takes to flush the carry-over gas out of the volume of one channel depends on the internal gas velocity:  Figure 5 shows that where the rotor enters the cold sector (for angles in the range 0 • to 5 • ), the temperature outlet profile shows a higher gradient. The reason for this is the carry-over gas from the hot sector. The time t f it takes to flush the carry-over gas out of the volume of one channel depends on the internal gas velocity: For the operating conditions in Figure 5, t f is about 10% of the time of a half cycle. Therefore, for almost the first 18 • during the rotation, just the carry-over gas from the process gas comes out to the supply air.
For numerical rotor simulation, some authors have used the approximated closed solutions method proposed by Kays and London [37] using the ε − NTU method, where the problem is simplified by assuming a constant heat transfer coefficient along the tube [13,14,17,20,21,23,[27][28][29]31,32,38]. Others have used finite difference methods for solving the same problem. In both cases, several authors have defined the axial heat conduction (in the direction of the flow) to zero, or radial heat conduction (perpendicular to the flow) to infinite. A detailed review on different approaches can be seen in the work of Klein and Eigenberger [30]. In this paper, it is assumed that thermal conductivity in the radial direction in the solid is infinite (which is the same to radial constant temperature). Additionally, the Biot number for the aluminum in the rotor conditions is in the order of magnitude of 10 −5 0.1, which is usually the limit for considering it as a lumped-system model Therefore the temperature gradients inside the aluminum are negligible.
Following the approach of previous literature [35,36], to solve Equations (5) and (6), these are discretized with central differences to a number of k + 1 grid points. Central differences have the benefit of being more accurate than forward or backward differences, since the truncation error is smaller. The time derivative is discretized according to the Euler implicit method, which is also known as the backward Euler method. This leads to a system of equations based on temperature differences in axial direction and in time. Common methods such as Gaussian elimination are constrained by numerical stability as the time step size ∆t decreases for a reasonable time resolution [39]. Therefore, the Newton-Raphson method is used. This method is also applicable when the physical properties, as in this case, depend on temperature. These equations were used in this work to develop an algorithm in C#. This method is robust in these cases and provides the aimed performance.
The design parameters have been incorporated for defining the step length ∆z, the hydraulic diameter d h , the void fraction ε, the volumetric specific surface of the rotor a v , and the inclination angle β. This data are required for the energy balance Equations (5) and (6), and the heat transfer coefficient models in Sections 2.1.1.4 and 2.1.1.5.
Different methods are used to calculate the different thermodynamic properties. Some data are based on correlations found in literature, others can be calculated using well-known thermodynamic laws or self-developed correlations based on literature data. A complete overview of the different methods is shown in Table 1. In the following section, we depict the different models used in this work for predicting the heat transfer coefficient α based on the Nusselt number Nu, and the advantages and disadvantages of each model are discussed. The Nusselt number is defined by: where the hydraulic diameter d h is defined by flow area A and perimeter P as follows: The above equation results in: To date, models for the heat transfer coefficient consider either constant heat flux transfer . q r to the wall in the radial direction (also called H boundary condition) or constant wall temperature ϑ S (also called the T boundary condition). The H boundary condition applies to counter-flow heat exchange, in which the wall heat flux to the fluid is constant along the z direction. The T boundary condition applies to constant temperature along the channel in the z direction but also periferically. These two boundary conditions serve as extreme cases. Actually, the rotor is expected to behave between these two boundary conditions.
Since the gas flows into the rotor in opposite directions in each sector of the housing, the temperature profiles along the channel inside the rotor are similar to those in a counterflow heat exchanger. However, this does not apply at the inlet and the outlet of the channel where the boundary conditions affect the temperature profile considerably, and in cases where the rotational speed n of the rotor is sufficiently slow, such that the rotor channels tend to reach the gas temperature. In these two cases, the heat transfer models should get closer to the T boundary condition.

Sine Duct Model
Shah and London [44] propose different values for Nu, depending on geometry.
In the H boundary condition, the authors also provide data for the H1 condition. H1 is a subcategory that applies to highly conductive materials, like copper and aluminum. The value of Nu H1 is always higher than the value of Nu T . The data from Shah and London [44] have the advantage that they refer to a similar geometry to that of the rotor. However, it only refers to fully developed flow.

Circular Tube Model
Gnielinski [33] reports two circular channels models: the constant heat flux model and the constant wall temperature model. The advantage of the circular channel model is that they are well studied and models for predicting the Nusselt number in the hydrodynamic developing entry region are available. Gnielinski [33] presents a model for constant heat flux, the H boundary condition, for circular tubes in developing flow, which allows calculation of the local Nusselt number.
This model is valid in the range 1000 > Pr > 0.7 and laminar flows. Since in this work the gas is air in the 273-350 K range, the Prandtl number is above 0.7. Since the internal gas velocity in the channel is below 5.5 m/s, and the diameter d h < 3 mm, the value of the Reynolds number is lower than 1000 at ρ G = 1.3 kg/m 3 and µ = 2.2 × 10 −5 Pa·s. The flow can be considered laminar due to Re < 2100 [27] in the aforementioned temperature range. Since both conditions are met, this model can be used.
In a similar way, Gnielinski [33] reports for constant wall temperature, the T boundary condition, which is valid for a wider range of the Prandtl number with Pr > 0.1.

Proposed Local Internal Nusselt Number for Sine Ducts (LINUS) Model
Since there is currently not a model in the literature considering both effects, i.e., geometry effect and hydrodynamic/thermal entry region effect, this study proposes a new model combining the advantages of both models reported above. The basic idea of this new model is to define the value of Nu on developed flow based on a correlation of the data gathered by Shah and London [44] that considers the geometry and to define Nu in the entry region based on the equations from Gnielinski [33] for each boundary condition. This proposed LINUS model is as follows: From the data from Shah and London [44], correlations using ς as defined in Equation (11) are obtained as shown in Equation (22) with the parameters from Table 2: In Figure 6 the Nusselt numbers obtained from these correlations are shown.
Energies 2020, 13, x FOR PEER REVIEW 10 of 26 Since there is currently not a model in the literature considering both effects, i.e., geometry effect and hydrodynamic/thermal entry region effect, this study proposes a new model combining the advantages of both models reported above. The basic idea of this new model is to define the value of on developed flow based on a correlation of the data gathered by Shah and London [44] that considers the geometry and to define in the entry region based on the equations from Gnielinski [33] for each boundary condition. This proposed LINUS model is as follows: From the data from Shah and London [44], correlations using as defined in Equation (11) are obtained as shown in Equation (22) with the parameters from Table 2: In Figure 6 the Nusselt numbers obtained from these correlations are shown. Based on the equations from Gnielinski [33], for the constant heat boundary condition is calculated as follows: where , in Equation (25) is calculated according to Equation (22) from the models from Shah and London [44]. Therefore, this model combines the Nusselt number for the sine geometry for developed flow from Shah and London [44] with the models for Nusselt number on the hydrodynamic developing region in Equations (23)-(25) from Gnielinski [33]. Based on the equations from Gnielinski [33], for the constant heat boundary condition Nu H1 is calculated as follows: Nu H1 = 3 Nu 1, H1 3 + 1 + (Nu 2, H1 − 1) 3 + Nu 3, H1 3 where Nu 1, H1 in Equation (25) is calculated according to Equation (22) from the models from Shah and London [44]. Therefore, this model combines the Nusselt number for the sine geometry for developed flow from Shah and London [44] with the models for Nusselt number on the hydrodynamic developing region in Equations (23)-(25) from Gnielinski [33]. Also based on the equations from Gnielinski [33], in the constant wall temperature boundary condition Nu T is calculated as follows: where Nu 1, T in Equation (28) is calculated according to Equation (22) from the models from Shah and London [44], following the same logic as Equation (25). Please note that in both boundary conditions, because it is based on the equations reported by Gnielinski [33], if z is high enough, the value of Nu tends to Nu 1 because Nu 2 and Nu 3 tend to zero. For example, in the operating conditions of this study if z > 55·d h , the effect of the entry region is negligible and the difference between the values of Nu and Nu 1 is less than 1%.
A comparison of the simulation results based on of these Nusselt number models (i.e., sine duct, circular duct or the proposed LINUS model with both boundary conditions) to experimental data is shown in the results section.

Pilot Plant Equipment
A pilot plant set-up was built for testing the simulation results in a laboratory at the Pforzheim University of Applied Sciences as shown in Figure 7. Filtered ambient air is blown into the rotor at low temperatures (from 0 to 25 • C) and then preheated by the rotor. Under normal conditions, the air would enter into an industrial process or a building for ventilation. In our case, Heater 2 simulates the heat that would come from an industrial process and provide hot air at temperatures up to 70 • C. This hot air is redirected into the rotor. The rotor is heated up by the hot air and the cooler air leaves the system. No purge sector is installed on the rotor. The whole system is shown in Figure 8 and a close-up to the heat exchanger at Figure 9. Also based on the equations from Gnielinski [33], in the constant wall temperature boundary condition is calculated as follows: where , in Equation (28) is calculated according to Equation (22) from the models from Shah and London [44], following the same logic as Equation (25).
Please note that in both boundary conditions, because it is based on the equations reported by Gnielinski [33], if is high enough, the value of tends to because and tend to zero. For example, in the operating conditions of this study if > 55 • , the effect of the entry region is negligible and the difference between the values of and is less than 1%. A comparison of the simulation results based on of these Nusselt number models (i.e., sine duct, circular duct or the proposed LINUS model with both boundary conditions) to experimental data is shown in the results section.

Pilot Plant Equipment
A pilot plant set-up was built for testing the simulation results in a laboratory at the Pforzheim University of Applied Sciences as shown in Figure 7. Filtered ambient air is blown into the rotor at low temperatures (from 0 to 25 °C) and then preheated by the rotor. Under normal conditions, the air would enter into an industrial process or a building for ventilation. In our case, Heater 2 simulates the heat that would come from an industrial process and provide hot air at temperatures up to 70 °C. This hot air is redirected into the rotor. The rotor is heated up by the hot air and the cooler air leaves the system. No purge sector is installed on the rotor. The whole system is shown in Figure 8 and a close-up to the heat exchanger at Figure 9.   Different devices are used in the pilot plant for measuring the following operating parameters, among them: • Temperature : 44 temperature sensors are distributed close before and after the rotor, with 11 on each side of the rotor in each sector. These groups of 11 sensors are evenly distributed along the rotating direction at a fixed radius from the center of the rotor, shown as crosses in Figure 10. Four additional temperature sensors are positioned some distance away from the rotor in each of the 4 sectors 1.1, 1.2, 2.1 and 2.2. These provide an average gas temperature. These 48 sensors are Pt-1000 sensors that use a 4-wire measurement method. These sensors are accurate between 0-150 °C as given by: In the temperature range of 0-70 °C, the maximum deviation of a single temperature sensor is less than 0.3 K. This allows the calculation of the thermal efficiency of the rotor from experimental results, at equal mass flows on both sides of the rotor, with a maximum deviation of ∆ = 0.8% from the temperature accuracy with . = 15 °C, . = 70 °C, and = 80%, which are the extreme conditions for the set-up.  Different devices are used in the pilot plant for measuring the following operating parameters, among them: • Temperature : 44 temperature sensors are distributed close before and after the rotor, with 11 on each side of the rotor in each sector. These groups of 11 sensors are evenly distributed along the rotating direction at a fixed radius from the center of the rotor, shown as crosses in Figure 10. Four additional temperature sensors are positioned some distance away from the rotor in each of the 4 sectors 1.1, 1.2, 2.1 and 2.2. These provide an average gas temperature. These 48 sensors are Pt-1000 sensors that use a 4-wire measurement method. These sensors are accurate between 0-150 °C as given by: In the temperature range of 0-70 °C, the maximum deviation of a single temperature sensor is less than 0.3 K. This allows the calculation of the thermal efficiency of the rotor from experimental results, at equal mass flows on both sides of the rotor, with a maximum deviation of ∆ = 0.8% from the temperature accuracy with . = 15 °C, . = 70 °C, and = 80%, which are the extreme conditions for the set-up. Different devices are used in the pilot plant for measuring the following operating parameters, among them:

•
Temperature ϑ: 44 temperature sensors are distributed close before and after the rotor, with 11 on each side of the rotor in each sector. These groups of 11 sensors are evenly distributed along the rotating direction at a fixed radius from the center of the rotor, shown as crosses in Figure 10. Four additional temperature sensors are positioned some distance away from the rotor in each of the 4 sectors 1.1, 1.2, 2.1 and 2.2. These provide an average gas temperature. These 48 sensors are Pt-1000 sensors that use a 4-wire measurement method. These sensors are accurate between 0-150 • C as given by: The pilot plant can work under the conditions shown in the following Table 3.

Research Design
In this section the experiments and simulations planned with the simulation software and pilot plant set-up as described in Section 2.1 are presented.

Validation Method
The thermal efficiency was tested in the pilot plant for a given empty tube gas velocity of ambient air . = 1.7 m/s and . = 4.2 m/s (which corresponds in mass-flux to 2.0 and 5.1 kg/(m 2 ·s) respectively) and different rotational speed between 4 and 24 rpm, which are the extreme values that can be reproduced at the pilot-plant.
The simulations were performed for the same operating conditions. The simulation used different models for the Nusselt number: Sine duct model [44], circular tube models [33], and the newly proposed LINUS model as described in Section 2.1.1.5 from this work. This procedure helps us to understand which models are more suitable for describing the experimental results.
Unfortunately, detailed data of the design and operational parameters of other rotors are rare in the literature, making it difficult to simulate rotors from literature with the developed simulation program. Therefore, the model that provides the best results compared to the experimental data from the pilot plant, is also used to simulate the data of rotor types published on the data from Eurovent [16], which have been the best example of detailed design information about rotors. The following rotor models in Table 4 are selected from the Eurovent database because they do not have a purge section that could potentially influence the results. These rotors are simulated under the Eurovent standard testing conditions of temperature, humidity, and empty tube gas velocity. In the temperature range of 0-70 • C, the maximum deviation of a single temperature sensor is less than 0.3 K. This allows the calculation of the thermal efficiency of the rotor from experimental results, at equal mass flows on both sides of the rotor, with a maximum deviation of ∆η = 0.8% from the temperature accuracy with ϑ 2.1 = 15 • C, ϑ 1.2 = 70 • C, and η = 80%, which are the extreme conditions for the set-up.

•
Volumetric flow . V: the volumetric flow was calculated based on gas velocity over a grid measurement according to DIN 16211 [44] using a hot-wire anemometer. With these values, the average flow velocity and the volumetric flow are calculated. With this method, the accuracy is around 9% of the measured value.
The pilot plant can work under the conditions shown in the following Table 3.

Research Design
In this section the experiments and simulations planned with the simulation software and pilot plant set-up as described in Section 2.1 are presented.

Validation Method
The thermal efficiency was tested in the pilot plant for a given empty tube gas velocity of ambient air v a2.1 = 1.7 m/s and v a2.1 = 4.2 m/s (which corresponds in mass-flux . m sp to 2.0 and 5.1 kg/(m 2 ·s) respectively) and different rotational speed n between 4 and 24 rpm, which are the extreme values that can be reproduced at the pilot-plant.
The simulations were performed for the same operating conditions. The simulation used different models for the Nusselt number: Sine duct model [44], circular tube models [33], and the newly proposed LINUS model as described in Section 2.1.1.5 from this work. This procedure helps us to understand which models are more suitable for describing the experimental results.
Unfortunately, detailed data of the design and operational parameters of other rotors are rare in the literature, making it difficult to simulate rotors from literature with the developed simulation program. Therefore, the model that provides the best results compared to the experimental data from the pilot plant, is also used to simulate the data of rotor types published on the data from Eurovent [16], which have been the best example of detailed design information about rotors. The following rotor models in Table 4 are selected from the Eurovent database because they do not have a purge section that could potentially influence the results. These rotors are simulated under the Eurovent standard testing conditions of temperature, humidity, and empty tube gas velocity. In the Eurovent data, the wavelength t * is defined slightly differently than in this paper: the wavelength is measured at the base of a wave on a curved arch inside the rotor [16]. The curvature, as explained already earlier in this paper, depends on the radius at which it is measured. This information is not specified by Eurovent data. For this reason, small differences are expected between the Eurovent results and the results of this paper, especially because the wavelength t * influences the wave height-to-length ratio and, therefore, also the Nusselt number.

Operating Parameters and Temperature Effectiveness
The operating parameters, which are empty tube gas velocity v a and rotational speed n, influence the thermal efficiency of the rotor. Therefore, the rotor efficiency was measured for a given empty tube gas velocity and different rotational speeds. This process was repeated at a different empty tube gas velocity.

Validation Results
The results of the simulations are the temperature profiles along the channel for the solid material and the gas for every time step on the rotating cycle as shown in Figure 11. In this figure, the hot gas is flowing from left to right in the channel and is heating up the aluminum matrix which is the solid phase. All the exhaust air temperatures provide the necessary information for calculating temperature effectiveness.  In the Eurovent data, the wavelength * is defined slightly differently than in this paper: the wavelength is measured at the base of a wave on a curved arch inside the rotor [16]. The curvature, as explained already earlier in this paper, depends on the radius at which it is measured. This information is not specified by Eurovent data. For this reason, small differences are expected between the Eurovent results and the results of this paper, especially because the wavelength * influences the wave height-to-length ratio and, therefore, also the Nusselt number.

Operating Parameters and Temperature Effectiveness
The operating parameters, which are empty tube gas velocity and rotational speed , influence the thermal efficiency of the rotor. Therefore, the rotor efficiency was measured for a given empty tube gas velocity and different rotational speeds. This process was repeated at a different empty tube gas velocity.

Validation Results
The results of the simulations are the temperature profiles along the channel for the solid material and the gas for every time step on the rotating cycle as shown in Figure 11. In this figure, the hot gas is flowing from left to right in the channel and is heating up the aluminum matrix which is the solid phase. All the exhaust air temperatures provide the necessary information for calculating temperature effectiveness.  Figure 11 shows that for the traditional rotor of 20 cm rotor length, from approximately 10% to more than 90% of the position inside the rotor channel, the temperature difference between the gas and the solid is approximately constant along the axial direc- Figure 11. Temperature profile along the rotor length (simulation result). Figure 11 shows that for the traditional rotor of 20 cm rotor length, from approximately 10% to more than 90% of the position inside the rotor channel, the temperature difference between the gas and the solid is approximately constant along the axial direction. Because the temperature difference is constant, the heat transfer in the radial direction from the gas to the wall is constant. Therefore, the H1 boundary condition is valid in this region. Boundary conditions govern the profiles of the heat transfer close to the inlet and outlet. In these locations inside the channel at 0-10% and 90-100% of the total length, the solid temperature is relatively constant. There, T boundary condition is valid, which has a lower value for Nu.
The pilot plant rotor was set to work under the following conditions shown in Table 5. The same parameters were used for a simulation using the different Nusselt number models. The results are presented in Figures 12-14. tion. Because the temperature difference is constant, the heat transfer in the radial direction from the gas to the wall is constant. Therefore, the H1 boundary condition is valid in this region. Boundary conditions govern the profiles of the heat transfer close to the inlet and outlet. In these locations inside the channel at 0-10% and 90-100% of the total length, the solid temperature is relatively constant. There, T boundary condition is valid, which has a lower value for . The pilot plant rotor was set to work under the following conditions shown in Table  5. The same parameters were used for a simulation using the different Nusselt number models. The results are presented in Figures 12-14.

Sine duct model from Shah and London
The following are the results from the test series on the pilot plant compared to the values predicted by the Shah and London model for heat transfer coefficient: The average difference between the experimental results and the H1 model of Shah et al. [44] is −2.0% points. In the pilot plant, increments of the rotational speed at the lower range produce a noticeable effect in temperature effectiveness, but increments at the higher range produce negligible increments in temperature effectiveness.  namic entrance length would be around 5 cm. This represents 25% of the duct length since the rotor length is 20 cm.

Circular tube model from Gnielinski
Here we present the results from test series on pilot plant compared to the values predicted by the models reported by Gnielinski [33] for the heat transfer coefficient for circular tubes. Figure 13. Results from pilot plant and simulation with circular tube models reported by Gnielinski [33].
The average difference between the experimental results and the T model from Gnielinski [33] is +3.8% points. This model shows the same slope of logarithmic growth as the model of Shah and London, but it shows much higher temperature effectiveness. Since an estimation shows that the heat losses in an extreme scenario would not be higher than 0.2% of the heat transferred in the rotor, the process can be considered adiabatic and it can be concluded that these models overestimate the heat transfer and do not match with empirical data. Furthermore, this can be explained by comparing the values of the Nusselt number for developed laminar flows reported by Gnielinski [33] ( = 3.66 and = 4.364) with the Nusselt number given by Shah and London ( = 1.95 and = 2.42). These results show also substantial differences between Nusselt numbers for sine ducts and circular tubes; thus the hydraulic diameter method is not applicable for rotor simulations with sine ducts.

 LINUS model
This model is based on the data from Shah and London for developed flow but it includes a prediction for the entry region based on the equations presented by Gnielinski Figure 13. Results from pilot plant and simulation with circular tube models reported by Gnielinski [33].
Energies 2020, 13, x FOR PEER REVIEW 17 of 26 [33]. The simulation results with the newly proposed LINUS model (with H1 and T boundary conditions, labeled hereafter as the LINUS model) are compared to the results from the pilot plant: The average deviation between the experimental results and the LINUS model with the H1 boundary condition is 0.1% points and the maximum absolute deviation was 1% point. This is a lower deviation compared to literature: the state-of-the-art simulation results had maximum absolute deviation of 3.5% points for Leong [12] or 3.71% on average for Özdemir [13]. This new model presents the same slope of logarithmic growth as the previous two models. The results from the simulation with the LINUS model match the pilot plant data. It is expected that the simulation results with the H1 and T boundary conditions are located above and under the pilot plant empirical data, respectively. This is because the temperature profiles (as shown in Figure 11) show that in some regions inside the channel, the H1 boundary condition is valid and in some others the T boundary condition is valid. However, the results are very close to the LINUS model with the H1 boundary condition; the reason for this is that the temperature difference between the solid and the gas is constant on almost 85% of the tube length (as shown in Figure 11). Since the temperature difference is constant, the constant heat transfer (H1 boundary condition) represents the reality better, as shown in Figure 14. Proposed Local Internal Nusselt

Sine duct model from Shah and London
The following are the results from the test series on the pilot plant compared to the values predicted by the Shah and London model for heat transfer coefficient: The average difference between the experimental results and the H1 model of Shah et al. [44] is −2.0% points. In the pilot plant, increments of the rotational speed at the lower range produce a noticeable effect in temperature effectiveness, but increments at the higher range produce negligible increments in temperature effectiveness. For both empty tube gas velocities (v a = 1.7 m/s and v a = 4.2 m/s), the resulting Nusselt numbers were Nu T = 1.95 and Nu H1 = 2.42. The calculation for the temperature effectiveness at the pilot plant shows an accuracy of approximately ±0.68%. For both boundary conditions, the model by Shah et al. [44] underestimate the heat transfer compared to the results observed in the pilot plant.
One possible reason could be that this model does not consider the effect of the heat transfer in the entry region. In this region, the heat transfer is higher than in the developed flow region and it, therefore, increases the average heat transfer coefficient throughout the entire channel length. This is especially true if the entry region length is significant compared to the overall channel length l. According to Shah et al. [44], the entrance length for laminar flow and sine ducts can be calculated as: For the rotor geometry installed in the pilot plant l hy + = 0.05. Since the Reynolds number in the pilot plant inside the channels can reach values up to 1000, the hydrodynamic entrance length l ent would be around 5 cm. This represents 25% of the duct length since the rotor length l is 20 cm.

Circular tube model from Gnielinski
Here we present the results from test series on pilot plant compared to the values predicted by the models reported by Gnielinski [33] for the heat transfer coefficient for circular tubes.
The average difference between the experimental results and the T model from Gnielinski [33] is +3.8% points. This model shows the same slope of logarithmic growth as the model of Shah and London, but it shows much higher temperature effectiveness. Since an estimation shows that the heat losses in an extreme scenario would not be higher than 0.2% of the heat transferred in the rotor, the process can be considered adiabatic and it can be concluded that these models overestimate the heat transfer and do not match with empirical data. Furthermore, this can be explained by comparing the values of the Nusselt number for developed laminar flows reported by Gnielinski [33] (Nu T = 3.66 and Nu H = 4.364) with the Nusselt number given by Shah and London (Nu T = 1.95 and Nu H1 = 2.42). These results show also substantial differences between Nusselt numbers for sine ducts and circular tubes; thus the hydraulic diameter method is not applicable for rotor simulations with sine ducts.

LINUS model
This model is based on the data from Shah and London for developed flow but it includes a prediction for the entry region based on the equations presented by Gnielinski [33]. The simulation results with the newly proposed LINUS model (with H1 and T boundary conditions, labeled hereafter as the LINUS model) are compared to the results from the pilot plant: The average deviation between the experimental results and the LINUS model with the H1 boundary condition is 0.1% points and the maximum absolute deviation was 1% point. This is a lower deviation compared to literature: the state-of-the-art simulation results had maximum absolute deviation of 3.5% points for Leong [12] or 3.71% on average for Özdemir [13]. This new model presents the same slope of logarithmic growth as the previous two models. The results from the simulation with the LINUS model match the pilot plant data. It is expected that the simulation results with the H1 and T boundary conditions are located above and under the pilot plant empirical data, respectively. This is because the temperature profiles (as shown in Figure 11) show that in some regions inside the channel, the H1 boundary condition is valid and in some others the T boundary condition is valid. However, the results are very close to the LINUS model with the H1 boundary condition; the reason for this is that the temperature difference between the solid and the gas is constant on almost 85% of the tube length (as shown in Figure 11). Since the temperature difference is constant, the constant heat transfer (H1 boundary condition) represents the reality better, as shown in Figure 14. Proposed Local Internal Nusselt Number for Sine Ducts (LINUS) model Given the accuracy of the equipment in this work, it is not possible to further improve the model and it is out of the scope of this project. For these reasons, the LINUS model with the H1 border condition is further used for simulations. The following Table 6 shows results of simulations of the LINUS model compared to the Eurovent data [16] of rotor models: These results show that the LINUS model is able to predict within 1.6% points the temperature effectiveness of these similar rotors published by Eurovent. The following aspects must be considered: (1) Eurovent published the data for rotor wave height with an accuracy of 0.1 mm. This means, for example, that wave heights from 1.55 to 1.64 mm, are reported as 1.6 mm. This is a major disadvantage for the analysis because according to the simulation, differences as small as 0.06 mm on wave height can already make a difference of 1.5% on temperature effectiveness. This could be a reason for differences between Eurovent empirical data and the simulation results of this work. A difference of 0.01 mm represents according to simulations a difference of approximately 0.1% points in temperature effectiveness. Therefore, wave height should be measured and reported with two decimal places (in millimeters) for achieving greater accuracy in the temperature effectiveness. (2) Eurovent empirical data have an accuracy of ±3% points [16]. The deviations seen in Table 6 between the simulation results with the LINUS model and the Eurovent data are significantly lower than the ±3% points of accuracy.
Nevertheless, the results are on average within 1.1% points of deviation from Eurovent empirical data. Given the accuracy of Eurovent data, there is good accordance of the simulation with experiments; therefore, these results validate the LINUS model for simulating rotors.

Rotational speed n
The experiments were performed under the conditions presented in Table 7. Since the H1 LINUS model showed the best accordance with pilot plant results, only this model was used.  Figure 15 shows that higher rotational speed provides higher temperature effectiveness. This effect is more significant at higher gas velocity. Higher rotational speed also has a drawback. There is a certain amount of gas that is enclosed within the rotor channels which is carried at the edge of each sector into the next sector (other side of the rotor housing).
If there is no purge-sector in the housing, ambient air (2.1) is carried over to the exhaust air (1.2) and process air (1.1) carried over to the supply air (2.2). This amount depends on the gas velocity v a and v in , the rotor length l, and the rotational speed n. The lower the gas velocity and the higher the rotational speed, the larger the amount carried over. This effect is outside the scope of this paper, but will be a possible future research aspect. According to pilot plant and simulation results, Figure 15 clearly indicates the decreasing improvement of the temperature effectiveness with increasing rotational speed. This effect is in agreement with results of other authors in Table 8. However, each author states a different number for this operational point. This might depend on the rotor geometry and empty tube gas velocity.

Empty Tube Gas Velocity v a
For the rotor of the pilot plant of this study (shown on Figure 9), an empty tube gas velocity of v a = 1 m/s (in mass-flux 1.2 kg·m −2 ·s −1 ) appears to be the optimum as in Figure 16. Other authors have also reported that a gas velocity close to 1 m/s provides better temperature effectiveness [18]. On the right side of the graph, it is clear that higher empty tube gas velocity v a decrease rotor thermal efficiency η in both simulation results and experimental data. Heidari-Kaydan et al. [18] and Sanaye [22] show the same tendency for higher gas velocity in their set-up. The reason is that with higher gas velocity the gas does not have enough time in contact with the surface for the heat transfer to take place. This tendency is shown also in the pilot plant results of this work on Figure 16. However, the pilot plant set-up used in this study cannot work on an empty tube gas velocity lower than v a = 1.7 m/s. Whether an empty tube gas velocity lower than 1 m/s decrease or increase the temperature effectiveness is a topic where there is no common agreement in previous Energies 2021, 14, 10 20 of 26 works: Heidari-Kaydan et al. [18] based on simulations report that empty tube gas velocity lower than 1 m/s would further improve the temperature effectiveness. Mahesh [26] by contrast claims that temperature effectiveness decreases in the range of low empty tube gas velocity (v a = 0.1 m/s). This same tendency is shown on Figure 16: lower empty tube gas velocity than v a < 1 m/s decrease the temperature effectiveness. The reason for this is that low empty tube gas velocity implies lower Reynold number and, therefore, lower Nusselt number in the entry region (see Equations (23)-(28)). Whether lower velocity than v a = 1 m/s increases or decreases temperature effectiveness is a topic that will be of interest for further research.
Energies 2020, 13, x FOR PEER REVIEW 20 of 26 velocity than = 1 m/s increases or decreases temperature effectiveness is a topic that will be of interest for further research. This simulated optimal empty tube gas velocity is = 1 m/s for the rotor configuration of this study. Since typical values for empty tube gas velocity in industrial applications are between = 2.7 and 4.1 m/s, a reduction of empty tube gas velocity would increase temperature effectiveness. To handle the same gas flow rate, the rotor diameter should be increased. An increase of rotor diameters to reduce empty tube gas velocity to = 1 m/s would increase the temperature effectiveness approximately by 10%.

Discussion
The vast majority of state-of-the-art simulations in this topic use either the ε-NTU method or computer fluid dynamics. This fact hinders the comparison of the results with other work. In this work, locally-calculated Nusselt number models with entry-region considerations were used. This is one of the main innovations of this work and what, according to the results, precisely increased the accuracy of the simulation models. Therefore, it is imperative that future work in this field report openly the deviations of the simulation to the empirical data. In the literature review, only few works reported the deviations between simulations and empirical data.
Additionally, constructive parameters and operational conditions should be transparently communicated, as in Tables 4 and 6. This is the absolute minimum to being able to simulate rotors. These are the essential parameters: rotor length , rotational speed , material thickness , wave height ℎ * , wave length * (or inclination angle * ), and empty tube gas velocity . These parameters have been shown to have a decisive role in the simulation results. Also the accuracy on the wave height ℎ * to ±0.01 mm proved to be necessary. Unfortunately, this is not disclosed with this accuracy in any literature.
An important aspect is that for the simulation homogeneous empty-tube gas velocity at the inlet sides of the rotor is assumed. In the pilot plant, this is true because the pilot plant is provided with flow straighteners. In reality, this is often not true and is an aspect  This simulated optimal empty tube gas velocity is v a = 1 m/s for the rotor configuration of this study. Since typical values for empty tube gas velocity in industrial applications are between v a = 2.7 and 4.1 m/s, a reduction of empty tube gas velocity would increase temperature effectiveness. To handle the same gas flow rate, the rotor diameter should be increased. An increase of rotor diameters to reduce empty tube gas velocity to v a = 1 m/s would increase the temperature effectiveness approximately by 10%.

Discussion
The vast majority of state-of-the-art simulations in this topic use either the ε-NTU method or computer fluid dynamics. This fact hinders the comparison of the results with other work. In this work, locally-calculated Nusselt number models with entryregion considerations were used. This is one of the main innovations of this work and what, according to the results, precisely increased the accuracy of the simulation models. Therefore, it is imperative that future work in this field report openly the deviations of the simulation to the empirical data. In the literature review, only few works reported the deviations between simulations and empirical data.
Additionally, constructive parameters and operational conditions should be transparently communicated, as in Tables 4 and 6. This is the absolute minimum to being able to simulate rotors. These are the essential parameters: rotor length l, rotational speed n, material thickness s, wave height h * , wave length t * (or inclination angle t * ), and empty tube gas velocity v a . These parameters have been shown to have a decisive role in the simulation results. Also the accuracy on the wave height h * to ±0.01 mm proved to be necessary. Unfortunately, this is not disclosed with this accuracy in any literature.
An important aspect is that for the simulation homogeneous empty-tube gas velocity at the inlet sides of the rotor is assumed. In the pilot plant, this is true because the pilot plant is provided with flow straighteners. In reality, this is often not true and is an aspect for further research on how inhomogeneous velocities could influence the temperature effectiveness. Additionally, the laboratory was leakage-free due to concrete measures and testing during the construction and installation of the system. This level of testing is probably not standard in simple heating, ventilation and air conditioning (HVAC) systems and leakages in heat recovery systems (HRS) could be more significant.
In this work, the average deviation of the simulation results (with the LINUS model and boundary condition H1) from the Eurovent empirical data is 1.1% points, and 0.1% points from the pilot plant empirical data. Hence, the simulation has been successfully validated, and the entry length effect played a key role. Previous work had higher deviations: a maximum absolute deviation of 3.5% points for Leong [12] or 3.71% in average for Özdemir [13] or were limited to test only one rotor exchanger [2,12,14,16,[19][20][21][22][23].
Rotational speed and empty tube gas velocity have an optimal operating point regarding temperature effectiveness. It will be further investigated whether an increase in temperature effectiveness actually reduces heating costs and CO 2 emissions. Furthermore, an optimization of the operational and design parameters shall be carried out in the future while considering energetic, economic and environmental aspects.

Conclusions
• Previous literature and state-of-the-art simulations of rotary heat exchangers had a limited accuracy: a maximum absolute deviation of 3.5% points for Leong [12] or 3.71% on average for Özdemir [13]. It was critical to consider the entry region effect on the gas phase for the simulation of rotors of 20 cm length and empty tube gas velocity from 1.7 to 4.2 m/s, since this was a key aspect to achieve higher accuracy of the simulation results.

•
Higher rotational speed up to 12 rpm provides higher temperature effectiveness.

•
Gas carry-over from the rotor should be further studied.
• Lower empty tube gas velocity (in the range from 1.7 to 4.2 m/s) increases the temperature effectiveness for the given rotor geometry. This was observed in the pilot plant and was predicted by the simulation. • Temperature effectiveness has a maximum at an empty tube gas velocity of v a = 1 m/s; whether lower velocity than v a = 1 m/s increases or decreases temperature effectiveness is a topic that will be of interest for further research.

•
Further studies shall determine the impact of increasing temperature effectiveness on manufacturing and operating costs and the impact on the environment.

•
An optimization of the operational and design parameters shall be carried out in the future while considering energetic, economic, and environmental aspects.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A. Thermophysical Properties
The following thermophysical properties from Table 1 are based on data from different authors. In this Appendix these data and the fitting curves are presented.

•
Aluminum density ρ S : This property was calculated according to the formula: where the K i values are given by:  And a comparison between the data from Rohsenow [40] (points) and the aforementioned equation (line) in the following Figure A1: Energies 2020, 13, x FOR PEER REVIEW 24 of 26 Figure A1. Aluminum density data from Rohsenow [40] and correlation (line). •

Aluminum thermal conductivity
This property was calculated according to the formula: where the values are given by: and a comparison between the data from Hatch [43] and the aforementioned equation in the following Figure A2: Figure A2. Aluminum thermal conductivity data from Hatch [43] and correlation. Figure A1. Aluminum density data from Rohsenow [40] and correlation (line).

•
Aluminum thermal conductivity λ S This property was calculated according to the formula: where the K i values are given by: and a comparison between the data from Hatch [43] and the aforementioned equation in the following Figure A2: Energies 2020, 13, x FOR PEER REVIEW 24 of 26 Figure A1. Aluminum density data from Rohsenow [40] and correlation (line). •

Aluminum thermal conductivity
This property was calculated according to the formula: where the values are given by: and a comparison between the data from Hatch [43] and the aforementioned equation in the following Figure A2: Figure A2. Aluminum thermal conductivity data from Hatch [43] and correlation. Figure A2. Aluminum thermal conductivity data from Hatch [43] and correlation.