Numerical evaluation of finite length tubes effects in Stirling engines heaters

The performance of Stirling engines is directly related to the amount of energy that, in the form of heat, participates in the thermodynamic cycle. A peculiar characteristic of this type of engine is the closed circuit of the working fluid, implying a periodic admission and extraction, ideally involving the whole working fluid, of the heat exchanged during the cycle without any exchange of material with the external environment. Several correlations have been proposed in the literature to predict the heat exchange in the hot side heat exchangers, but most of them are based on nondimensional numbers, usually expressed in terms of an oscillatory Reynolds number and a non-dimensional length, that even if they take into account the amplitude and frequency of the oscillating flow inside the tube with respect to its diameter, do not include any dependency upon the length of the tube. Nevertheless, the length of the tube can have a great impact on the performance of a Stirling engine. The heater forms a significant part of the dead volume of the engine, making the optimization of the volume to surface ratio necessary. The friction losses increase by increasing the length of the tubes, determining a negative impact while the exchange surface increases. Even the Nusselt number in the inner side of the tubes changes along the length, achieving the largest values alternatively in the first portions of the tube lengths because of entrance effects (Graetz problem). The picture is made even more complex because of the special velocity profiles that develop in oscillating flows in tubes. One of the effects less investigated in previous studies, which we could call a “breathing effect”, regards the amount of working fluid that, in a real engine, can effectively travel from the hot side to the cold side, thus reaching the conditions of nominal heat exchange with the thermal source and sink of the cycle. An idealized configuration has been devised to investigate, using CFD simulations, these effects. Results reporting how the friction coefficient and the Nusselt number depend on the finite length of the tube will be illustrated. * Corresponding author: francescosaverio.marra@stems.cnr.it © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). E3S Web of Conferences 313, 03003 (2021) https://doi.org/10.1051/e3sconf/202131303003 19° International Stirling Engine Conference


Introduction
Most models for the prediction of the performances of Stirling engines rely on the availability of empirical correlation for the estimate of the performances of single critical components: regenerator and heat exchangers are probably the best examples [1]. Unfortunately, a large set of available data in the specific operating conditions of Stirling engines, firstly the reciprocating regime, is not available. In the past years, the major research effort was devoted to the regenerator, due to the large influence of this component on the attainable performances when a proper design is conducted [2]. Less attention was devoted to the heat exchanger, partly because this aspect seemed better understood in the enormous body of literature devoted to this physical process [3], partly because its influence and the peculiar behavior in the not so common oscillating (or reciprocating) regime [4] was widely recognized only later. It is only in recent years that papers addressing this issue have appeared in the literature.
Several researchers investigated heat transfer under an oscillatory flow regime. One of the first thorough investigations of the heat transfer phenomena in the tube of the heater of a Stirling engine was conducted by [4]. They pointed out that for several Stirling engine analyzed, the working fluid does not pass entirely through heat exchangers during a cycle, making the inclusion of the effect of finite length of the tube necessary for a correct prediction of the heat flux. Kornhauser and Smith [5] recognized several phenomena due to the complex interaction of flow displacement, temperature, and pressure that become not synchronized when subject to the driving force of the pistons moving with a phase shift. The importance of the effect of the phase angle between the angular velocities of the pistons was confirmed by Kanzaka and Iwabuchi [6,7] which determined an optimal phase shift for the heat transfer of 180 deg. Shin and Nishio [8] analyzed the flow and heat transfer in a pipe with heated and cooled short sections traversed by an oscillating fluid.
Among the several results, they pointed out that for amplitudes of oscillations smaller than the heating/cooling lengths, the averaged heat transfer coefficients are close to those for the unidirectional, fully developed laminar flow. An interesting study on the effect of several parameters (including the tube length) on the heat transfer was performed by Worlikar and Knio [9], focused on the cooler of a thermos-acoustic refrigerator. They measured the critical length at which no further addition of heat transfer occurs by increasing the length of the tube beyond the critical one, correlated to the amplitude of fluid displacement. They also recognized the importance of the flow structure at the end portion of the tubes. The complex field of the main variables influencing the heat transfer, velocity, temperature, and pressure, in the heater tubes was experimentally investigated also by Bouvier et al. [10]. They found a large difference between the local and the globally averaged Nusselt numbers, meaning that a single section evaluation, as done in many experimental setups for the measurement of the heat transfer, cannot be reliable. Kuosa et al. [11] performed an extensive comparison of correlations for both unidirectional and oscillating flow heat transfer and proposed some measures for heat transfer enhancement. To measure the heat transfer, they adopted the mechanical structure of a real Stirling engine, thus including the effects of the phase shift of the movement of the pistons driving the oscillating flow. Xiao et al. [12] obtained the velocity profile in the fully developed oscillating regime for values of the dynamic Reynolds number spanning both laminar and turbulent regimes. A fixed length test section was considered and no compressibility effects were taken into account. In a companion paper [13], other effects were included, principally the distribution of the heat flux among the several tubes that form a real heat exchanger. They highlighted that the unidirectional correlations for the Nusselt number are inadequate to predict the performances of heat transfer in oscillating flows, especially at high Reynolds numbers. Ni et al. [14] focused on the real shape of the tubes of the hot side heat exchanger, which often are U-shaped instead that purely linear. Based on measurements performed on an apparatus based on a single piston to drive the flow oscillation, thus without compression effects, they proposed different correlations for helium, nitrogen, and carbon dioxide as working fluids. Barreno et al. [15] developed new time dependent correlations for the Nusselt number in presence of oscillating flow, based on numerical simulation results with a constant property fluid. Their findings indicate an increase in the Nusselt number due to the reciprocating flow. Ilory et al. [16] used a single acoustic driver to promote the oscillatory flow and obtain details of the velocity and pressure profiles in the test section. The focus was on the high frequency of the oscillation for the application to thermo-acoustic engines. Xin et al. [17,18] investigated the enhanced heat transfer in the cooler pipes of Stirling engines, showing the possibility to increase the heat flux by introducing a helical wire along the internal wall of the heat exchanger tubes. However, the interaction with the flow developing for the driving action of piston and displacer was not taken into account.
Not many papers addressed the question of whether an optimal configuration of the oscillating flow in the tubes of the heat exchanger exists. The first is the contribution of Zarinchang and Yarmahmoudi [19], which determined theoretically the optimal values of the number of tubes, the length, and the diameter for an alpha type Stirling engine. Yin and Ma [20] derived an analytical solution for an ideal incompressible configuration, identifying the importance of the Prandtl number. Lombardi et al. [21] using an optimization procedure, based on empirical correlations, modified the original design of an alpha-type Stirling engine heater subject to high thermal flux in a fluidized bed, obtaining equivalent performances reducing the number and length of the tubes. However, all these results are obviously affected by the validity of the empirical correlation adopted. Patil and Gawali [22] gave one of the best contributions to this question. They, studying the heat transfer in the oscillating flow of water, recognized the existence of an optimal frequency depending on the tidal displacement length. Recently, Islas et al. [23] performed an extensive optimization study in which length, diameter, and number of the heater tubes were taken into account. Their study indicated that diameter and number of heater tubes are among the parameters that mostly influence the efficiency of the engine.
Despite the previously recalled contributions that appeared on this subject, several effects depending on the specific conditions induced by the Stirling thermodynamic cycle have been not examined. Indeed, the correlation parameters adopted, the kinetic Reynolds number Re and the dimensionless oscillation amplitude A [24], refer to the simpler circumstance of an oscillating flow far from the entrance length of the pipe, and with synchronous pistons movement, so that compression and expansion effects are not included. It is worth noting that the first hypothesis mentioned leads also to the exclusion of the effect of the heat exchanger pipe length and the relative volume of the heat exchanger (part of the dead volume in the Stirling engine nomenclature) to the volume displaced by the moving parts. Further effects depend on the shape of the inlet/outlet sections of the heat exchanger pipes, the non-uniform exposition to the external heat flux along the pipe surface, the effect of the conjugate heat transfer in the solid wall of the pipe. These effects are closely related to the specific engine configuration and manufacturing and coupling design with the heat source, and therefore they are not considered in this work.
This contribution aims to highlight some of the aspects not completely included in the previous studies. The focus is on the hot side heat exchanger, assuming that on the cold side, thanks to a favorable interaction with liquid water, proper dimensioning is less critical. In addition, the analysis is limited to pipe heat exchangers, one of the most widely adopted shapes of heat exchangers in Stirling engines. The aim of this work is also to contribute to a better understanding of the several parameters influencing the performances of heat exchangers of Stirling engines in a real environment, by varying the portion of the tube exposed to the heat flux. Then, the effects of the tube length and the tube portion exposed to the external heat flux at varying the Reynolds number and the oscillation amplitude are investigated by means of CFD (Computational Fluid Dynamics) numerical simulations. After a brief excursus of the characteristic relevant parameters of existing Stirling engines, the adopted model is introduced, the selection of the conditions explored in the analysis is presented, and the results are then discussed. Some conclusions are eventually draft from the most relevant results.

Reference parameters
In order to generalize the solutions of specific cases, it is usual to introduce reference quantities that produce non-dimensional groups whose values identify similar solutions. The approach proposed by Zhao and Cheng [24] has been adopted and it is here briefly recalled to clarify the relationship between the non-dimensional groups and the dimensional parameters adopted on the simulations.
A first reference quantity is the average distance that a fluid element travels in an infinite pipe of diameter during a full oscillation. Please note that in a Stirling engine, this distance is usually larger than the length of the pipe of the heat exchanger: this means that the fluid participating in the heat exchange can be not uniform, depending on the mixing occurring in the piston's chambers. This length, normalized by the pipe diameter, defines the first non-dimensional group A : A second reference quantity is introduced after the velocity of the oscillation. Let indicates the angular velocity of the oscillation [rad/s], a Reynolds number can be defined as A third parameter is defined combining the two previous ones: This parameter helps to identify the regime of motion in the pipe [24,25]: for ≤ 350 the flow in the pipe can be assumed laminar, while for ≥ 3000 the flow can be assumed fully turbulent. Two different regimes are identified in the range 350 < < 3000: for 350 < < 750 the flow is not anymore fully laminar, but the formulation adopted is still the laminar one; for 750 < < 3000, where the flow becomes transitional, the turbulent formulation is adopted.
A fourth dimensionless parameter is introduced to describe the effect of the penetration length of the fluid in the tube of the heat exchanger. This is computed as the ratio of the oscillation amplitude on the tube length : It represents the number of sweeps of the tube by the semi-oscillation of the fluid. Table 1 reports some of the geometrical and operating parameters, relevant for the current study, of Stirling engines collected from the literature. Even if these data represent only a small number of exemplar real engines, it is interesting to note that the values of Re and A looks clearly correlated: despite the possibility to select independently each of them, practical design considerations seems to impose a constraint requiring an almost linear dependency, at least for the lower values of A , that remain valid for the parameter while seems to assume a plateau for the parameter Re , as recognizable in Figure 1.

. Governing equations
The compressible formulation of the Navier-Stokes equations is adopted for the numerical modeling of the phenomena. They include the balance equations for total mass, momentum ad energy. For the variable volume domain Ω( ), indicated with the velocity of points on the boundary ( ), they can be written as: Momentum Balance In these relations, is the density of the working gas, the velocity of a fluid element, = , with representing the pressure, the identity tensor, and representing the stress tensor for a Newtonian flow with dynamic viscosity , is the gravity acceleration, the total energy and the heat flux. For our problem, the driving force is the pressure exerted by the forced movement of the pistons, so that gravity acceleration is neglected. Also neglected are heat fluxes apart from the Fourier contribute, so that = ∇ , where is the temperature. Coupling these equations with the ideal gas state equation = and the caloric state equation: we obtain the closed set of governing equations.
In the case the flow cannot be assumed laminar, a turbulent flow formulation is adopted. In this case, all state variables appearing in Eqs. 6-8 are assumed Favre time averages of instantaneous variables and the phenomenological coefficients and are corrected by adding the contributes and of the turbulent fluctuations. These contributions are computed from the knowledge of the turbulent variables (turbulent kinetic energy) and (dissipation rate of turbulent kinetic energy) that comes by adding the respective balance equations [35]. The so-called Scale Adaptive Shear Stress Transport formulation is here adopted [36], due to its better performance with unsteady flow configurations in presence of adverse pressure gradient. Indeed, adverse pressure gradients develop at each inversion of the flow in the oscillating regime. The computational domain is illustrated in Figure 2. The problem formulation is completed by assigning the initial and boundary conditions. Fluid at rest in the position of 0 crank angle (displacer in the top dead center) at the cold temperature and uniform charging pressure is selected as the initial condition to mimic the practical start-up stage of the Stirling engines. Boundary conditions correspond to no-slip conditions on the whole boundary, different fixed hot temperature/mixed adiabatic conditions on the connecting pipe between the two cylinders, a fixed cold temperature on the ring between the pipe and the lateral walls of the "piston" cylinder, and adiabatic conditions on the "displacer" cylinder walls. These conditions will be specified later. Finally, a sinusoidal movement, with a phase angle of 90 degrees, of the bottom boundaries of the two cylinders is assigned.

Numerical Solution
The solution to the problem is obtained using the open-source code OpenFOAM [37]. This suite of libraries, applications, and solvers, allows building the discretization mesh, defining the problem in the discretized formulation, and obtaining the numerical solution. Details of this procedure are not given here, being behind the scope of this paper and because they are described in the user guide of the software. It is here reported that the grid adopted is fully 3-dimensional, with size in the and coordinates (see Figure 1) for the fixed section of 110 and 10 cells respectively. Further cells are adopted for the moving part of the cylinders, depending on their required extension to keep a similar accuracy. The radial mesh size is 15 cells for the pipe section and 30 cells for the cylinders respectively. Being the results described thereafter just preliminary ones for this investigation, a complete grid sensitivity analysis has been not yet completed. However, by comparison of the numerical solutions with the experimental results reported in [24], this grid size appeared adequate to get a good approximation of the unsteady solution, at least for the velocity field, see Figure 2. Other important assumptions are relative to the form adopted to discretize the energy equation, that is rewritten adopting enthalpy as state variables, and the expressions adopted to compute the phenomenological coefficients. A constant value of is assumed (1007 kJ/kg K). The dynamic viscosity is computed as a function of temperature with the Sutherland correlation while the thermal conductivity is computed using the Eucken approximation [38].

Numerical setup
The realization of the different conditions to explore is realized adopting the simple domain geometry reported in Figure 2. Three interconnected cylinders form it. The two at the extrema, named D (Displacer) and P (Piston) and both of radius , have the first section from the closed end of length variable with time, to drive the oscillating motion, and the section adjacent to connecting pipe fixed to account for a dead volume and modify the compression ratio. The central cylinder, of reduced radius , works as a connecting element but also as a heat receiver. Heat rejection occurs on the circular crown of cylinder P.
A synchronized stroke of the end face of both cylinders determines, in the case of incompressible flow, a bulk movement of the gas inside the connecting pipe of amplitude given by = / . This is the value adopted in the definition of the nondimensional group A . In order to keep constant the heat rejection surface, the dimensions and are fixed. Therefore, to change the parameter A it is varied the length of the stroke, following the relation: The selected Reynolds number is assigned, fixed the working gas, the initial average pressure (thus defining a reference value for the kinematic viscosity = / and the pipe diameter, by changing the value of the rotation speed : Air is adopted as the working fluid at an initial average pressure of 1 bar. ( 1 3 ) where is the tube diameter and is the average thermal conductivity of the gas in the section , evaluated, in conformity with the values used in the CFD computations, as a function of temperature following the Eucken approximation adopting the values ( , ). Finally, a global average Nusselt number is defined for the entire pipe averaging the contributions during the cycle of the oscillations:

Definition and evaluation of the Nusselt number
Several correlations are available in the literature for the estimate of this important parameter. Barreno et al. [15] have reviewed several past correlations and proposed an advanced correlation that allows computing the instantaneous global Nusselt taking into account, in addition to A and Re , also the ratio / . In the following, the averaged values in a cycle computed adopting this correlation will be indicated as Nu .

Results
The convergence of the numerical simulation to a fully developed flow required a long computation time: many cycles are required to reach the equilibrium between the heat addition along the tube and the heat subtraction on the annular wall of the cooler cylinder. The development of such conditions was ensured by the analysis of time signals of temperature and pressure at several locations in the computational domain. Figure 3 shows the time profiles for temperature and pressure in the case with Re = 302, A = 21, and = 1.4. About 1 s of physical time is required to get, using these parameter values, almost fully developed, periodic, conditions. However, all simulations have been continued for a significant time to ensure the full development of the periodic regime, and the averaging procedure is performed only on the last period computed.  Figure 4 shows the velocity and temperature fields along the tube for the same computation. 8 equally time-spaced subsequent intervals in a cycle period are reported. Some interesting features emerge. During the first half period, the flow occurs mainly from the hot side on the left (displacer) to the cold side on the right (piston). During this stage, the penetration of the hot gases, together with the rise of the temperature due to compression, reduces the temperature difference with the wall, progressively limiting the heat transfer. Inversion of the flow movement occurs firstly along the tube walls, causing the well-known annular flow of the oscillating regime.  This is somewhat different from the flow evolution observed in the case of synchronous movement of the driving pistons, and thus in absence of compression, as recognizable in Figure 5. A steep traveling front in the temperature of the gases establishes which pushes, in both directions, the gases ahead (see Figure 5 left). This effect impacts the distribution of the heat flux, as shown in Figure 5 right, that especially in the forward phase, is shielded by the temperature distribution that realizes along the radius, reported in Figure 6. The lowest Nusselt numbers are observed at the instants 2/8 , when almost all the tube is filled with the hottest gases and the axial velocity shows an almost parabolic profile, and 7/8 , when again the axial velocity reassumes a parabolic profile for the entire length of the tube and the normal derivative of the temperature at the wall is the minimum.

Effect of
An increase of A , keeping constant Re and the tube length, equals to modify the number of sweeps of the flow in the tube (realize deeper breathing at the constant frequency). This effect has been investigated adopting the parameters reported in Table 2. A linked effect is to increase the parameter and thus shift towards the turbulent regime. In the current preliminary results of this investigation, the transition to turbulence has been promoted at the highest values of A by doubling the parameter Re to avoid an excessive deformation of the moving section of the mesh that exhibited numerical instabilities. The transition to the turbulent regime appears to occur for 695 < < 869, in accordance with the previous findings of Akhavan et al. [25] and Zhao and Cheng [24]. . Red circles, Nusselt computed with the correlation by Gnielinski (see Barrenoet al. [15]). Blue plus, Nusselt computed with the correlation by Barreno et al. [15]. Re = 302 for A < 80. Re = 604 for A > 80.
The values of the Nusselt number computed in this work are significantly lower than those predicted by the correlation proposed by Gnielinsky [15] and even than those predicted by Barreno et al. [15], as reported in Figure 7, a discrepancy that needs further investigation. It is important to note that both correlations were derived by assuming imposed temperatures at both sides of the tubes and constant thermal fluid properties. In the present work, no boundary conditions were imposed at the entrance and exit of the pipe, being both sections part of the computational domain. This difference has important consequences: firstly, as already described by Simon and Seume [4], during the cycle a two-zone heat transfer mechanism occurs. A first zone at contact with the wall, where heat transfer can be hindered, during part of the cycle, by the hot annular flow of gases coming from the hot side. A second zone, the inner bulk that can also flow in a reverse direction with respect to the annular external region, thus limiting the heat transfer along the main flow direction. A further effect superposes to the previous one at higher values of : a compression wave travels along the tube causing the maximum temperature at the wall to be locally larger than the wall temperature, in conjunction with the effect due to the compression of not-in-phase pistons movement, as documented by the results reported in Figure 8 for Re = 302, A = 40 and = 2.7.

Effect of but at constant .
In order to get a constant number of sweeps while changing the parameter A it is necessary to change the length of the tube as indicated in Table 3. The corresponding computed Nusselts are plotted in Figure 9. These values are different from those reported in Table 2 for the same values of A and Re , indicating the need to take into account both the length of the tube and the number of sweeps.

Effect of the portion of the tube exposed to heat flux
The last comparison is done between the results already discussed, all obtained by exposing all the surface of the pipe to a constant hot temperature, and the results obtained by exposing to a hot temperature only 1/3 of the entire length located on the side exposed to cooling. For the remaining length, the pipe surface is assumed adiabatic. A graphical comparison with the corresponding values for the full length of the pipe heated is reported in Figure 9. In both cases, the Nusselt number is computed assuming all the surfaces of the pipe as participating in the heat exchange with the hot source of the Stirling engine. This assumption corresponds to the practical situation in which only a portion of the heater is effectively exposed to the nominal heat flux, due to installation requirements, inhomogeneous heat source, etc. Obviously, the so computed Nusselt numbers are much lower. However, a beneficial effect is observed in which the reduction is not proportional to the portion of tube length exposed. Indeed, as reported in the last column of Table 4, the ratio is always greater than the ratio of the surface extension exposed to heat flux that is 0.33. This result is explained by the circumstance that the entrance region of the pipe contributes with much higher values of the Nusselt number, an effect that is present also in the section separating the heated from the adiabatic portion of the tube, building a thermal entrance effect, as clearly shown in Figure 10 for the case Re = 302, A = 21 and = 1.4 . Fig. 10. Computed and Nu along the pipe at 8 equally spaced instants during a cycle. Re = 302, A = 21, and = 1.4 but only 1/3 of the tube surface exposed to heat flux.

Conclusions
A new configuration is proposed to study the many effects that several geometric and operational parameters can have on the effective heat transfer in the hot side heat exchanger of a Stirling engine formed by a bundle of tubes. The preliminary results, here presented, illustrate the effect of the finite length of the tube, its relation with the limited fluid displacement, and the coupling with the mixing occurring in the other volumes that forms the overall space occupied by the working fluid.
Two main conclusions can be draft from the obtained results. The first is a prediction of lower value of averaged Nusselt numbers than those reported in the literature; a result that needs further investigation and possibly an experimental confirmation. The second is the possibility to design heat exchangers with significantly shorter tubes than those usually adopted, especially if the heat transfer mechanism on the external side of the tube is enhanced by means of careful coupling with the heat source.
This study is still far to be complete. Apart from the necessary validation of the results here reported against experimental data, other effects need to be included, as a different portion of the heated part of the tube, for instance centrally placed or close to the full adiabatic volume of the hot side, and a different length of the fixed sections of the displacer and piston volumes, which correspond to a different total dead volume and compression ratios.
Luigi Acampora acknowledge partial financial support with CNR grant on Research Project DIT.AD017.012 -Processi e Tecnologie per l'Energia da Fonti Rinnovabili. Computational resources have been gratefully provided by University of Sannio, Project G.E.M.M.E., Prof. Gaetano Continillo.