A 3D-CFD study of a γ-type Stirling engine

Accepted Manuscript A 3D-CFD study of a γ-type Stirling engine Lukasz Kuban, Jakub Stempka, Artur Tyliszczak PII:

S0360-5442(18)32375-2

DOI:

https://doi.org/10.1016/j.energy.2018.12.009

Reference:

EGY 14262

To appear in:

Energy

Received Date: 22 August 2018 Revised Date:

28 November 2018

Accepted Date: 2 December 2018

Please cite this article as: Kuban L, Stempka J, Tyliszczak A, A 3D-CFD study of a γ-type Stirling engine, Energy (2019), doi: https://doi.org/10.1016/j.energy.2018.12.009. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

A 3D-CFD study of a γ-type Stirling engine Lukasz Kuban1 , Jakub Stempka, Artur Tyliszczak

RI PT

Czestochowa University of Technology, Faculty of Mechanical Engineering and Computer Science, Al. Armii Krajowej 21, 42-201 Czestochowa, Poland

Abstract

M AN U

SC

The goal of this work is to perform detailed 3D CFD analysis of heat transfer and flow dynamics inside a commercial gamma-type Stirling engine, including its structural details and variability of the working spaces. The study focuses on an impact of different initial engine parameters, i.e., filling pressures, rotational speeds and heater temperatures, on the engine’s characteristics. One of the biggest advantages of using CFD in such kind of problems, is that it allows to reduce models uncertainty resulting from inaccurate assumptions inherent in the low-order models (analytical formulas for heat transfer coefficient or pressure losses etc.). Such empirical formulas are based on the literature data, therefore they are valid only for some range of engines’ operational conditions. The results were validated against available experimental data, showing good agreement in terms of overall trends. Comparison of different turbulence and regenerator models demonstrated their little impact on the engine indicated work. The applied CFD model allowed us to accurately and deeply analyse 3D unsteady flow features including the velocity field, vorticity, temperature, friction coefficient and production/dissipation of the turbulent kinetic energy. We focused on identification of recirculation/stagnation regions as well as the regions of increased turbulence intensity.

TE D

Keywords: Stirling engine, heat exchange, CFD/turbulence modelling, dynamic meshes

1. Introduction

AC C

EP

Protection of the natural environment and preservation of the worldwide energy resources require efficient use of all kinds of fuels, including fossil fuels, bio-fuels as well as green fuels (wind, water, sun). These issues are usually connected with an impact of energy production technologies on the surrounding world, and particularly, with combustion processes and unavoidable problems of toxic gas emission and pollution. Awareness of irreversible degradation of the environment caused by quickly growing industry in 60s-70s led to environmental law regulations. According to them industrial furnaces and gas turbines installed in large and medium size power plants had to be redesigned and optimised. This had been happening over the last 20-30 years when large, single chamber burners evolved into compact, smaller multi-burner configurations [1] using new ideas for low emission combus1 Email

address: [email protected] (L. Kuban)

Preprint submitted to Energy

tion technologies. These include approaches based on reheat systems [1, 2, 3], a colourless distributed combustion (CDC) [4], MILD combustion [5, 6, 7], an enhanced mixing [8, 9, 10, 11], reverse [12, 13, 6] and forward flows configurations [14, 15]. These works show a huge interest and progress that have been made approaching eco-friendly combustion technologies. The achievements are undisputed, but unfortunately, it seems that in this field there is not much left to do in the future, unless some breakthrough will occur. Modern prototypes of the combustion chambers/burners dedicated to oxycombustion systems [16, 17], and particularly, their future combination with hydrogen fuel leading to a perfectly clean combustion with the water vapour as the only combustion product [18], seem to be promising, but they are still under development. The aforementioned technologies have one common feature, i.e., at the end of the process there is some energy (heat) left, which is transferred back to the environment, for instance, as a warm water from power plants or with exhaust gases from December 3, 2018

ACCEPTED MANUSCRIPT

erly analysed using these models. Nevertheless, those models are used for preliminary engine design and optimization purposes [29, 30, 31, 28, 32]. It is generally agreed, that modelling of complex heat and mass transfer processes such as occurring in the Stirling engines require much more sophisticated tools. With the advent of powerful computers the research on the Stirling engines can be performed using CFD tools. By far the CFD is the only approach, which besides being capable of modelling conjugated flow and heat processes, also allows to take into account real geometries of the engines, including their moving parts. In contrary to the analytical methods mentioned previously, the CFD tools can handle complex processes that are often neglected or simplified in loworder models (e.g., laminar-turbulent transition). Although, CFD allows detailed study of the local flow characteristics it also requires some modelling assumptions to be made. Nevertheless, those assumptions are certainly less restrictive, easily perceivable and rather oriented on specification of the flow type (laminar vs. turbulent, and if so, which turbulence model), boundary conditions (e.g. constant or spatially/temporarily varying wall temperature), medium properties (e.g. constant or temperature dependent viscosity and thermal conductivity). From the numerical point of view, since the flow scales occurring the Stirling engine vary significantly, its subsequent cycles are unsteady and not fully repeatable, therefore, modelling is very challenging task. It gets even more complex when moving elements such as power piston and displacer are taken into account, requiring use of time-varying computational domains and moving meshes. Moreover, the CFD results may be sensitive to the choice of numerical schemes, discretization and integration methods, time step size and computational mesh quality, so improper selection of any of them may lead to stability issues or lack of convergence. The successful modelling of the configurations like Stirling engines using CFD methods requires not only the deep understanding of heat and mass transfer processes but also significant experience in numerical simulations in order to obtain meaningful results free from misinterpretation. Although the CFD results presented in literature are often obtained from the computations performed for significantly simplified cases, they enable much more advanced analyses than the analytical models. For example, two dimensional (2D) simulations of a β-type engine have been performed

AC C

EP

TE D

M AN U

SC

RI PT

chimneys. There are many ways aiming at retrieving wasted thermal energy, e.g., regenerative heat exchangers, Rankine cycle-based power plants, devices using thermoelectric effect and the Striling cycle-based thermal engines, which are getting more attention from the industry and research centres over the last few decades. There are three basic configurations of the Stirling cycle-based engines, the so-called α-type, β-type and γ-type that differ in the placement of underlying working volumes or piston types [19]. All variants of the Stirling engines share similar advantages, the most important are good efficiency, smooth operation and easy maintenance [19, 20, 21]. The one specific advantage, that makes them exceptional among other heat engines, is that the Stirling engine can operate on any kinds of heat sources, not only those originating from the combustion processes. For example, the so-called low-temperature-differential (LTD) Stirling engines [22, 23, 24] are regarded as one of the most promising solutions to the global warming problems. They can work in the range of temperature differences below 100◦ C therefore, they can be combined with solar, biological and geothermal heat sources. The main disadvantage, which preclude the use of LTD Stirling engines on large scale is related to the fact that an increase in the power output requires larger pressure what makes sealing more difficult. Additionally, an issue relating all the Stirling engines types is that their power output relay solely on the conditions of an external heat source, so the engines are basically slowly responsive and hard to control. Design and optimization of the Stirling engines are always proceeded by analytical and numerical modelling using tools with growing level of complexity and accuracy. In that sense, the models are classified as of the first-, second- and thirdorder. The analytical models are practically zeroand one-dimensional in which the entire engine is divided into axially connected sections (usually 3 or 5). Among the first order models an improved isothermal analysis formulated by Gustav Schmidt in 1871 is the one most commonly used [25, 26, 27]. Second and third order models take into account transient heat transfer [28] and thermal an mechanical losses [29] but at most they can result in the variations of thermal and flow variables only in the axial direction. Neither variations in the transverse direction nor effects caused by geometrical details (displacer dimensions, compression/expansion chamber size, dead zone volume, etc.) can be prop-

2

ACCEPTED MANUSCRIPT

ling engines, even if applied with significant geometrical simplifications and assumption of the flow laminarity. An important observation is drawn from the work, i.e., the temperature, velocity and heat transfer show non-uniform 3D distributions, what suggests that: (i) the complete engine has to be considered, rather than its quarter; (ii) as many as possible geometrical details should be included as they may significantly affect the flow field. Recently, Alfarawi et al. [37], using COMSOL Multiphysics software [36] performed the numerical simulation of a γ-type high differential temperature Stirling engine (ST05-CNC) (similar one as in the current study). In the simulations the flow was assumed to be laminar and the moving engine parts and changing volumes of the cylinders were handled using Lagrangian-Eulerian (ALE) approach. Although, the authors used simplified 2D domain, it included almost every geometrical detail of the engine. This work was claimed as a first attempt to simulate of such complex Stirling engine configuration using CFD tools. The achieved results highlighted an importance of optimisation of the phase angle and the size of the dead spaces. An advanced 3D CFD parametric study of an α-type Stirling engine was recently performed by Almajri et al. [40]. They analysed the effects of various design parameters and operating conditions (regenerator porosity, charge pressure, matrix wire diameter and hot/cold ends temperatures) on the engine performance in terms of indicated power output. Cheng and Chen [41] carried out 3D simulations of a prototype 1 kW β-type Stirling engine and based on a measured engine power output, they performed validation of the numerical model. Similarly as Almajri et al. [40], Cheng and Chen also studied impact of different values of working parameters (filling pressure, heating temperature, regenerator porosity, etc.) on the engine P-V diagrams. To the best knowledge of the authors, the only work which is devoted to the fully 3D simulation of a Striling engine (particularly a γ-type engine) was done by K. Mahkamov [38]. It concerned the engine working in combination with an industrial biomass gasifier. A sophisticated algorithm was applied to deal with the changing engine volumes and the number of cells in the computational meshes varied during the engine cycles. The remeshing was achieved by adding or removing the cell layers above the displacer and power piston during the strokes. The used meshes were relatively coarse, as could be seen in the figures presented in [39]. Neverthe-

AC C

EP

TE D

M AN U

SC

RI PT

by Salazar and Chen [33] with the use of their own numerical code based on a multi-block approach and structured meshes. The flow was assumed laminar and axisymmetric. During every cycle meshes were modified by expanding or compressing their volumes. Although, this approach does not seem to be the best option, as the local mesh resolution, and thus a solution accuracy, had been changing in time, depending on the position of piston and displacer, they managed to accurately predict an exemplary solution. The authors analysed the temperature and velocity distribution as well as the time variations of the heat flux along the walls. They observed highly non-uniform temperature contours and occurrence of the recirculation regions. It was found that impingement of the flow on the engine walls is responsible for the major heat transfer mechanism. This research was further extended to studies of low-temperature-differential γ-type engine focusing on characteristics and intensity of the heat transfer in a very simple threedimensional (3D) configuration [34]. The computational model assumed geometrical symmetry, and hence, only a quarter of the real engine was analysed. Again the laminar flow model was used and the meshes were modified in time in the same manner as for the above discussed 2D case. It was confirmed that in the 3D model the heat transfer is the strongest also in the region of the flow impingement. The results revealed temperature and velocity variations in the azimuthal direction leading to occurrence of three dimensional vortical structures visualised by the velocity contours. The same engine was analysed in [35, 24] where alteration of selected geometric and operational parameters were investigated for an impact on the engine performance. The research included variability of the engine speed, temperature difference, phase angle between power piston and displacer (φ), gap between displacer and displacer cylinder (δ) and also length of linkage bar (ld ). Among others, it was found that by reducing δ by nearly half, the indicated power increased by 32%. They observed that with increasing ld , the power and efficiency increased linearly by 10.08% and 13.83%, respectively. Concerning the variations of φ only small changes have been reported, however, it was shown that φ = 100◦ yields the maximum indicated power, whereas φ = 90◦ guarantees the best efficiency and therefore this value was advised as the preferred one. The authors demonstrated a huge potential of the CFD tools for analysis and optimization of Stir-

3

(a)

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

Figure 1: CAD geometry of the ST05G Stirling engine with the names of its main parts (a) and sample results showing temperature distribution on the walls of the heater, regenerator and cooler (b).

AC C

EP

TE D

less, the obtained velocity field allowed identification of recirculation and stagnation zones. Results revealed that the velocity field differs significantly from one tube to another, despite the same geometry. Such differences may lead to a strongly nonuniform thermal load on the heater parts and therefore, its inefficient usage. The above review shows that the numerical studies of the Stirling engines are advancing but a vast of analysis are performed for simplified configurations. The complex multidimensional and transient simulations are essential for understanding and optimisation of inner engine processes. A rapid increase of computing power and maturity of CFD tools allow performing computational simulations using more physical models. The present paper is aimed at the application of CFD tools to perform a 3D-time dependent Unsteady Reynolds Averaged Navier-Stokes (URANS) simulation with moving and deforming meshes of the commercially available Stirling engine (model ST05G) [42] designed by the German constructor Dieter Viebach. This engine was analysed experimentally in many research for verification and improvements of the analytical models [28, 29]. The goal of this work is twofold: (i) to perform URANS simulations of the ST05G us-

4

ing different turbulence models and regenerator parameters without introducing significant geometrical simplifications; (ii) to show complexity, multidimensionality and non-linearity of the heat and mass transfer processes occurring in its working spaces. We focus on 3D unsteady flow features and identifications of the recirculation/stagnation regions as well as the regions of increased turbulence intensity. Beside that, we analyse the impact of the rotational speed, pressure load and heater temperature on the engine power. To our best knowledge the results of simulations of a practical Stirling engine have never been presented in such details as in this work. The paper is organised as follow, the next section introduces the description of analysed engine and its operating parameters. The mathematical model, the numerical method and applied meshes are presented in section 3. The obtained results, are firstly compared with data from reviewed literature, next 2D/3D plots presenting flow characteristics are given and described in Section 4, which is followed by conclusions.

ACCEPTED MANUSCRIPT

2. ST05G Stirling engine characteristics

Table 1: Main engine parameters.

Operational parameter

value

unit

The unit studied presently is a gamma type Stirling engine ST05G of Viebach design (Fig. 1a). Its main parameters are briefly summarised in Table 1. The dimensions of the engine are approximately 410×375×780mm and its total weight is equal to 33.5kg. The engine consists of two cylinders which axes are coplanar and perpendicular. On the top of one of the cylinders (named hereafter ”hot cylinder”) a heater consisted of 20 stainless steel 8×1mm tubes are mounted. Just below the heater around the cylinder, the regenerator assembly is placed. The regenerator is made of steel wire mesh of diameter 0.05-0.06mm and porosity 0.9, its length is equal to 53.5mm. Further below the cooler part is placed. The cooler itself consists of two main parts, i.e., the one which is visible comprises an external aluminium wall while the inside inset is made of a machined piece of aluminium cylinder to the form of a finned heat exchanger which external and internal diameters are equal to 124 and 101mm, respectively, and its length is equal to 84mm. Both cylinders are mounted into cuboidal aluminium crankcase. The power piston is placed in the other cylinder (named hereafter ”cold cylinder”) on top of which an aluminium head is attached that is connected with the “hot cylinder through a stainless steel tube. The tube creates a cooled duct for working gas and conists of two pipes, i.e., the visible one (35×1.5mm) which is just an external wall and inside one (28×1.2mm) which confines the flowing gas. In between, the flow passage for cooling medium is designed. The second end of the tube is mounted below the cooler and the flow is directed towards space under the displacer. The role of the displacer is to exert only a small amount of pressure on the gas providing its flow through the heater, regenerator and cooler, this piston does not produce any work. For minimizing the thermal load on displacer walls, near the heater, displacer has 156mm long hollow head which diameter is equal to 98mm that is made of thin (0.8 mm thickness) stainless steel. The second piston, the one that produces work, which diameter is equal to 84.5mm and length is equal to 76.8mm, is made of aluminium alloy. The stroke of both pistons is constrained by the crank mechanism and is equal to 75mm. The flywheel is made of cast alloy and comprises a mechanical energy reservoir. The engine nominal rotational speed is equal

power output (at 500RPM)

500

[W]

max. speed

750

[RPM]

max. heater temperature

550

[◦ C]

working gas pressure

10

[bar]

stroke length

75

[mm]

piston bore diameter

84.5

[mm]

displacer bore diameter

98

[mm]

SC

RI PT

2.1. Unit description

M AN U

to 500RPM at which it produces 500W, assuring the right temperature of heat sources, for which the higher one is limited to 550◦ C. Such a performance is obtained for the nitrogen as a working gas, however, other gasses are also possible to use.

AC C

EP

TE D

2.2. Working principle of the Stirling engine Every Stirling engine works based on a transformation of the thermal energy into the mechanical one. Figure 1b shows a sample distribution of the wall temperatures of the main parts in the engine. The movement of the power piston is caused by the gas expansion/compression occurring when the gas is heated/cooled flowing through the heater/cooler. The power piston is situated within ”cold cylinder” and the displacer is placed in the ”hot cylinder”. The phase shift between the piston and displacer is constrained by the geometry of the crank mechanism. It is equal to 90 degrees what corresponds to 1/4th of the working cycle. The crank mechanism converts linear motion of the pistons to rotations of the shaft. The engine kinematics, including exact positions and velocities of the moving parts are fully determined in time by two variables, i.e., shaft angle and rotational frequency. Every working cycle of the engine consists of a few characteristic points which may be described as follow: a) when the displacer reaches the ”bottom dead” point, the power piston is on the half-way to it; after passing the ”bottom dead” point the displacer starts to move back; b) when the power piston reaches the ”bottom dead” point, the displacer is on the half-way to the ”top dead” point; after passing the ”bottom dead” point the power piston starts to move back;

5

ACCEPTED MANUSCRIPT

Table 2: Analysed cases.

RPM[min−1 ]

T[◦ C]

P[bar]

R600 P10 T773

500

10

600

R600 P10 T673

400

10

600

R600 P10 T573

300

10

600

500

6

600

500

3

600

500

10

360

500

10

180

R600 P6 T773 R600 P3 T773 R360 P10 T773

SC

R180 P10 T773

RI PT

Case name

compression; afterwards the whole cycle starts again.

M AN U

Figure 2: Schematic representation of the Stirling engine cycle.

3. Modelling approach

c) when the displacer reaches the ”top dead” point, the power piston is on the half-way to it; afterwards the displacer starts to move down; d) when the power piston reaches the ”top dead” point and the displacer is located in a half-way of the cylinder.

TE D

After another quarter of the cycle, the process starts over.

EP

From the thermodynamic point of view a practical Stirling engine ST05G does not follow the theoretical Stirling cycle but rather its adiabatic version (see Fig. 2), which could be described as follows: 1-2) pseudo-isochoric heat addition; the compressed working fluid flows through the regenerator and absorbs heat on the way to the expansion space;

3.1. Computational domain and working conditions Figure 3 presents the computational domain, which consists of cylinder volume, the heater tubes, the tube connecting both cylinders and the regenerator (modelled as a porous medium, which parameters are discussed later). The velocities of the displacer and power piston are determined by the value of a rotational speed of the engine. Knowing that the phase shift between the power piston and displacer is equal to φ = 90◦ they velocities are defined as: Vdisp = ω · Lrod · sin(ω · t)

(1)

Vpist = ω · Lrod · sin(ω · t − π/2)

(2)

AC C

where ω = RPM × 2 · π/60 is angular velocity, t is the time and Lrod = 0.0375m is the half of the displacer and power piston swept length. We analyse three different temperatures of the heater (cooler temperature is set to 37◦ C for all cases), three initial pressure loadings and three rotational speeds as shown in Table 2. All heat sinks and heat sources were considered as isothermal and all the remaining wall boundaries were assumed to be adiabatic. Therefore, the heat losses to the environment were neglected in the present simulations. The air with constant heat capacity cp = 1006J/kg/K and with the molecular viscosity µ(T ) changing according to the Sutherland law was used as a working fluid. The thermal conductivity was computed as κ(T ) = cp µ(T )/Pr with the Prandtl number equal to Pr = 0.7.

2-3) pseudo-isothermal expansion; the expansion space is heated externally, and the gas undergoes expansion; the total engine’s volume increases and the pressure decreases; 3-4) pseudo-isochoric heat removal; the gas passes through the regenerator, hence some of its heat is transferred to the regenerator for the use in the next cycle; 4-1) pseudo-isothermal compression; the working medium flowing to the compression space is cooled down in the cooler and then undergoes 6

ACCEPTED MANUSCRIPT

(3)

ef ∂τij ∂ ρ¯u ¯i u ¯j ∂ p¯ ∂ ρ¯u ¯i + =− + ∂t ∂xj ∂xi ∂xj

(4)

SC

¯i ∂ ρ¯ ∂ ρ¯u + =0 ∂t ∂xi

RI PT

3.2. Mathematical model The flow inside the engine was assumed to be compressible, turbulent and non-isothermal. We applied the URANS approach with the realizable k −  turbulence model to compute the Reynolds stresses [45]. The continuity equation, the NavierStokes equations and the energy equation in the framework of URANS are given as:

M AN U

  ¯ ¯ + p)¯ ∂ ρ¯E ∂(¯ ρE uj ∂ ∂ T¯ ef + = (κ + κt ) + τij u ¯i ∂t ∂xj ∂xj ∂xj (5) where the effective viscous stress tensor is defined as:    ∂u ¯i ∂u ¯j 2 ∂u ¯l ef τij = (µ + µt ) + − δij (6) ∂xj ∂xi 3 ∂xl In the above equations the bar symbol denotes the Reynolds time-averaging procedure, while the remaining symbols denote: t - time, xi - spatial direction, ρ - density, ui - velocity component, p pressure, E - total energy, µ - molecular viscosity, κ - thermal conductivity. Equations (3)-(5) are complemented with the equation of state for perfect gas p¯ = ρ¯RT¯, where R is the individual gas constant. The variable µt , is turbulent viscosity, which needs to be modelled and κt = Cp µt /Prt is the turbulent thermal conductivity, where Prt = 0.85 is the turbulent Prandtl number. The standard version of the k −  model was formulated by Launder and Spalding [43] in 70s. It comprises of two equations, one for the turbulent kinetic energy k and the second for its dissipation rate , which are related to the turbulent viscosity as: k2 (7) µt = ρCµ  where Cµ in the standard k −  model is treated as a constant value equals 0.09. Since its invention the k −  model and their later modifications ( k −  realizable and RNG) were extensively used in many practical applications and they are probably ones of the most used models within RANS approach. In the k −  realizable model used in the present work the Cµ is no longer constant, it is computed

AC C

EP

TE D

Figure 3: Computational domain of the ST05G Stirling engine.

depending on the mean strain rate, rotation rate, the angular velocity of the system and local values of k and  [45]. The assumption for derivation of the k −  family models is that the flow is fully turbulent. In theory, these models should not be applied for laminar or transient regimes, however, practice shows that these limitations are often neglected as it is difficult to assess whether the flow in whole domain is turbulent or laminar. The transport equations for k and  are given as [44, 45]:

∂ ∂ ∂ µt ∂k (ρk) + (ρk¯ ui ) = [(µ + ) ]+ ∂t ∂xi ∂xj σk ∂xj (8) + Gk + Gb − ρ − YM

∂ ∂ ∂ µt ∂ (ρ) + (ρ¯ ui ) = [(µ + ) ]+ ∂t ∂xi ∂xj σ ∂xj + ρC1 S − ρC2 7

2  √ + C1 (C3 Gb ) k k + ν

(9)

ACCEPTED MANUSCRIPT

where  C1 = max 0.43,



TE D

M AN U

SC

RI PT

p η k , η = S , S = 2Sij Sij η+5  (10) and Sij is the strain rate tensor. Additional terms in Eq. (8) and (9) denote: Gk - generation of turbulence kinetic energy due to mean velocity gradients, Gb - generation of turbulence kinetic energy due to buoyancy, YM - contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. The definitions of these terms can be found in [44]. The symbols σk , σ , C1 , C2 , C3 are the model constants, which were calibrated based on various flow types including shear flows, homogeneous flows and decaying isotropic turbulence. Their typical values are [44, 45]: C1 = 1.44, C2 = 1.9, σk = 1.0 and σ = 1.2. Compared to the standard k −  model its realizable version has superior ability to capture the complex flow structures involving separation, rotation and even recirculation regions. It can be also used to model boundary layers under strong adverse pressure gradients. As will be shown later all these phenomena occur in the inner volumes of the Stirling engine and therefore it can be assumed that the chosen turbulence model is appropriate for the present studies. However, in order to assess the level of influence of the turbulence modelling on the results additional test computations using k−ω SST model have been performed for a few selected cases. Details of this model can be found in [46, 44].

the hexahedral cells were used wherever it was possible without significant deterioration of the mesh quality (see Fig. 4). All remaining parts, which geometrical complexity did not allowed to use structural meshes, were built using unstructured tetrahedral elements. For instance such places are the connections between major parts presented in Fig. 5. The mesh in the ”hot cylinder” was built using three sub-zones, this simplification allowed to reduce complexity of cell addition and removal during the piston motion. The separate sub-zones were created for heater pipes and the regenerator. The cooler was subdivided into three sub-zones. To sum up the whole mesh was built from 8 structural regions of which 5 were dynamically modified and 2 unstructured regions. The computations were performed assuming the transient regimes with the time step equal to ∆t = 0.2 × 10−3 s. The number of time steps per single cycle were equal to 500, 833 and 1666, depending on the rotational speed and cycle period respectively. In between the successive steps, due to the motion of piston and displacer the mesh size varies dynamically from 9.2 to 11.6 million cells. The spatial discretisation is based on the 2nd order upwind scheme for the Navier-Stokes and energy equations and the 1st order upwind scheme for the turbulent kinetic energy and its dissipation. The time integration was performed applying the 1st order implicit Euler approach, which in the ANSYS Fluent v.17 is the only option for the cases with moving elements and deforming meshes. Within each time step the scaled residual dropped to 10−4 in less than 20 subiterations for each equation. The convergence rate was only slightly dependent on the working phase of the engine and direction of the movement of piston and displacer. The computations have been performed on CPU Cluster using 24 Intel Xeon E5-2680 cores. After an initial transient stage dominated by strongly unsteady effects occurring when the engine start (e.g. temperature or pressure jumps) the computations of every cycle took 15, 22 and 43 hours respectively for engine speed equal to 600, 360 and 180RPM. The solutions in terms of velocity field, pressure temperature and some derived quantities (e.g. vorticity, heat fluxes or friction coefficients) were monitored as averaged quantities over the cycles, either in the whole volume of the domain and at selected locations. The regenerator was modelled as a porous body with the porosity equal to α = 0.9 with vis-

3.3. Computational mesh and solution algorithm

AC C

EP

The set of equations (3-9) is solved using the ANSYS Fluent v.17 software [44]. Unlike in the research performed in [33, 24], where the number of cells in the meshes were kept constant and they dimensions varied in time, in the present work we could benefit from the use of the commercial software and ready-to-use sophisticated algorithms for remeshing the volumes around the moving pistons. When expansion or compression space volume increases, the leading cell layer expands. When the cell height increases beyond value of 0.5mm, the whole layer gets split and new layer which continues to expand is being created. Analogically, when the volumes of working spaces decrease, given mesh layer gets squeezed and eventually merges with neighbouring layer. In order to be able to use the feature of dynamic mesh layering the mesh had to be carefully built. Structural meshes based on 8

RI PT

ACCEPTED MANUSCRIPT

SC

Figure 4: View of the computational mesh in the displacer cylinder in the region close to the bottom wall.

AC C

EP

TE D

M AN U

cous and internal resistances set respectively to 9.433e+09 m−2 and 76090m−1 (taken from [47]). Effects of the regenerator porosity were studied by Cheng and Chen [41] who analysed its impact on engine output power in the range α = 0.7 − 0.95. They observed only 3% change of the engine power output (1800W - 1850W) with the maximum at α = 0.9. Hence, knowing that the porosity has only minor impact on the results and exact parameters of the resistance were not given for the modelled unit, in the present studies we have analysed an impact of the viscous and internal resistances of the porous material. The analysed test cases included decrease and increase of regenerator viscous and internal resistance by 20%.

Figure 5: View of the computational mesh below the displacer in the region close to the bottom wall and between the power piston and the cooler pipe.

4. Results Firstly, in this section we concentrate on the accuracy of performed simulations and compare the obtained results with available experimental data from literature. Unfortunately, the reference data were not obtained for exactly the same engine unit or working conditions as used in the present studies. In many experimental studies the engines were slightly modified from one experiment to another. In such cases the comparisons have mainly qualitative character and we focus on similarity of the observed trends, e.g. increase of the engine power vs. increasing heater temperature or filling pressure. 4.1. Numerical model validation The engine was initially filled with the air at the temperature equals to 320K, the initial pressure was equal to pload = 10bar and the temperature of the

9

ACCEPTED MANUSCRIPT

4.1.1. Pressure and volume variations in one cycle Figure 6 shows the evolutions of the volume and pressure in the function of the crank angle. The total volume refers to the complete volume of the engine (working and dead spaces). The expansion and compression space are presented separately and corresponds to the ”hot” and ”cold” cylinders. The pressure lines represent the mass weighted average values of the absolute pressure in the whole engine. Analysis of the volume changes reveal the phase shift between the expansion and compression space volumes. The maximum pressure occurs at the crank angle equal to 210 degrees. At this particular position the compression space reaches its minimum volume. One can also notice the difference in the maximum pressure in the compression and expansion spaces. This difference is partially due to pressure losses in the regenerator, as discussed by Hachem et al. [29], but also because the gas in the compression space has lost some heat in the cooler, and hence its pressure is lower.

Work [J]

Power [W]

R600 P10 T773

66.5 (67.1)

665 (671)

R600 P10 T673

51.1

511

R600 P10 T573

32.2

322

R600 P6 T773

42.3

423

R600 P3 T773

23.6

236

R360 P10 T773

64.2 (63.3)

385.2 (379.8)

R180 P10 T773

59.9 (57.9)

179.9 (173.7)

M AN U

Total Exp. Space Comp. Space

14

3

13

1.2E-03

12

6.0E-04

Pressure [bar]

1.8E-03

Volume [m ]

SC

Case

RI PT

Table 3: Calculated indicated work and power for the studied cases. The values in the brackets refer to the results obtained with k − ω SST model.

0

60

120

180

TE D

11

240

300

360

Crank Angle [deg]

EP

Figure 6: Volume and pressure changes for whole engine and for expansion and compression spaces as a function of crank angle for the case R180 P10 T773 .

AC C

heater was Th = 773K (see Table 2). While the engine was running the mean temperature and the mean pressure increased. After the initial transient phase, which took approximately 20 cycles, the working conditions of the engine stabilised and every subsequent cycle looked virtually the same, except for small deviation resulting from an unpredictable nature of turbulence. We noted that achieving different working parameters could be obtained either by starting from the initial conditions and following the above procedure or by changing the working parameters when the engine is running. We tested both methods and it turned out that the final (mean) working cycle was independent on the way in which this state was reached.

4.1.2. Cycle-to-cycle variations As already mentioned the heat and fluid flow inside the engine are strongly unsteady and some differences between subsequent cycles can be expected, for instance, caused by turbulence. Sample results illustrating such a behaviour are presented in Fig. 7 showing temporal variations of the average temperature monitored between the expansion space and the inlet/outlet regenerator surfaces. Figure 7 shows also as the average heat flux on the heater and cooler walls. These results were obtained for three cases (R180 P10 T773 , R360 P10 T773 , R600 P10 T773 ) using two turbulence models (realizable k −  and k − ω SST). Careful analysis of the solutions shows that the temperature and heat flux changes in subsequent cycles were not identical. This phenomenon can be observed even after long simulation time and this confirms a non-negligible impact of turbulence on unsteadiness of the successive cycles. For this reason in the analyses related to the global engine characteristics, e.g., to the data concerning P-V plots or engine power the representative results were obtained by averaging the instantaneous solutions over several periods. 4.1.3. Impact of turbulence modelling and regenerator parameters Analysis of the results presented in Fig. 7 shows that the turbulence modelling approach has only

10

RI PT

ACCEPTED MANUSCRIPT

(b) R360 P10 T773

M AN U

(c) R600 P10 T773

SC

(a) R180 P10 T773

Figure 7: Temperature and heat flux variations in the selected engine parts for the cases: R180 P10 T773 , R360 P10 T773 , R600 P10 T773 computed using the k −  realizable model (lines with ) and k − ω SST model (lines without symbols).

TE D

minor impact on the obtained instantaneous solutions. One can observe that the differences are small, only quantitative and decreases with the increasing RPM speed since the flow velocity and turbulence intensity increases. Further comparisons will be presented on the P-V plots analysed in the following subsections.

AC C

EP

Impact of varying regenerator parameters discussed in Section 3.3 is presented in Fig. 8 that shows the temporal variation of the temperature and heat flux in the selected engine parts for the reference case R600 P10 T773 computed using the k −  realizable model. As can be seen in the presented plots in comparison to the differences resulting from the application of various turbulence models, the impact of decreasing or increasing the porous material resistance seems to be larger. This is clearly visible if one considers the temperature variations in the heater expansion space and the heat flux on the heater walls. Compared to the reference case the increase of the resistance of the porous material causes the periodic rise of the temperature. For small time span, approximately 0.04s, it increases even above the heater wall temperature (Th = 773K), which results in the negative heat flux on the walls. The heat flux in the heater for higher resistance is most of the time lower than in the reference case. On the contrary, for the decreased resistance the temperature in the heater is

Figure 8: Temperature and heat flux variations in the selected engine parts for the reference case R600 P10 T773 computed using the k− realizable model for varying regenerator parameters: solid lines - reference case; lines with 4 - porous material resistance increased by 20%; lines with  - porous material resistance decreased by 20%.

smaller compared to the reference case and the heat flux is always positive. On the other hand, the heat flux on the cooler walls is apparently unaffected by the regenerator parameters. As will be shown later these seemingly evident instantaneous differences have very small influence on the global engine characteristics. 4.2. Experimental verification The results in terms of the indicated work and power output obtained for all cases analysed in this work are given in Table 3. The numbers in the brackets correspond to the solutions achieved with 11

SC

RI PT

ACCEPTED MANUSCRIPT

(b)

M AN U

(a)

Figure 9: PV-diagrams for: a) the same heater temperature, different rotational speeds and turbulence models; b) the k −  realizable model , the same rotational speed and different heater temperatures.

AC C

EP

TE D

the use of the k − ω SST model and the remaining ones are for k −  realizable. It can be seen that for both cases the results strongly depend on the working conditions. However, the differences between the solutions obtained applying k − ω SST and k− realizable model are very small and as previously observed (Section 4.1.3) they decrease with increasing rotational speed. Regarding the effective work the relative difference between the models calculated as |Workk− − Workk−ω |/Workk− × 100% for 180RPM is equal to 3.4% and it drops down to less than 1% for 600RPM. Hence, it seems that the turbulence model has a small impact not only on the instantaneous solutions, as could be seen in Section 4.1.3 but also on the global characteristics of the engine. In the following parts, we will analyse their dependence on the heater temperature, the initial pressure load and the rotational speed of the engine.

4.2.1. P-V diagrams The P-V diagrams obtained for different working and modelling parameters are presented in Figs. 9 and 10. Note, that the values of the work presented in Table 3 are computed as the surface integrals of the P-V plots. The results shown in Fig. 9 are presented separately for the cases with: a) different rotational speeds (180, 360 and 600 RPM), two turbulence models and the same heater tempera-

Figure 10: PV-diagrams for the reference case R600 P10 T773 with varying regenerator parameters.

12

ACCEPTED MANUSCRIPT

ture Th = 773K; b) different heater temperatures (Th = 573K, 673K and 773K), k −  realizable model and the same rotational speed (600RPM). Note that in all the cases (a and b) the initial pressure was the same pload =10bar. It can be seen in Fig. 9a that for different engine speeds the diagrams are similar, but not the same, as one could expect. Interesting is that increasing the rotational speed the higher values of the maximum pressure and also larger work per cycle (see Table 3) are reached. Note that the same tendency is observed when the k − ω SST model is applied. Compared to the solutions obtained using the realizable k −  model the plots are slightly shifted up, which means that the mean pressure is higher. However, as could be seen in Table 3 the resulting indicated works were similar. The reason for such somewhat surprising dependence of the solutions on the engine speeds can be related to the fact that for larger values of rotational spped the flow velocity, and thus the Reynolds number (Re), increases. During the post-processing of the results, the heat transfer is computed based on the effective viscosity including both the physical and turbulent one. The latter increases when the turbulence is stronger, which usually happens for larger flow velocity. Hence, the observed shift in the P - V diagrams can be attributed to over-prediction of the turbulent viscosity by both turbulence models. Certainly, this aspect has to be verified in future analyses. Impact of the heater temperature on the P-V diagrams is shown in Fig. 9b. It is clearly seen that providing more thermal energy to the engine translates to the larger area in between P-V lines, and thus the larger useful work. For the heater temperatures 773K, 673K and 573K the calculated work per cycle is equal to 66.5J, 51.1J and 32.2J, respectively. An increase of the temperature by 17% (from 573K to 673K) leads to 58% increase of the indicated power, then the increase of the temperature from 673K to 773K (15%) results in another 30% increase of the indicated power. Dependence of the P-V diagrams on the regenerator parameters is shown in Fig. 10 which presents the result corresponding to the solutions previously shown in Fig. 8. There, the differences in the temporal behaviour of temperature and heat fluxes could be easily noticed, at least for some engine regions and for small periods of time. Apparently, as can be seen in Fig. 10 these differences have practically no impact on the P-V plots, neither on the resulting effective work. When the porous material

AC C

EP

TE D

M AN U

SC

RI PT

resistance is decreased the effective work equals to 66.8J, whereas for the increasing resistance it drops down to 65.9J. Compared to the reference case for which the work is equal to 66.5J (see Table 3) the relative difference is smaller than 1%. As discussed before, such a relatively small impact of the regenerator parameters were previously reported by Cheng and Chen [41], where by changing the regenerator porosity by 20% the maximum differences in the results reached at most 3%. Similar conclusions can be drawn from [49], where the effectiveness of regenerators built from various materials (aluminium, stainless steel, copper and monell) were compared experimentally for the ST05G engine. It was found that the regenerators built from the cooper and stainless steel characterized visibly higher efficiency. However, for these cases for given operating conditions (pload =8bar and Th =773K) the difference in the engine power was smaller than 6%. This was somewhat surprising taking into account that the differences between the physical properties of cooper and stainless steel are as high as 20% in thermal capacity, 14% in density and 15 times thermal conductivity. Hence, taking into account the present solutions and literature data one can assume that the applied modelling approach and parameters of the regenerator should not introduce considerable uncertainty into the results.

13

4.2.2. Parametric analyses Impact of the heater temperature on the indicated power (P ) is presented in Fig. 11 along with the results of experimental work of Alfarawi et al. [37], Hachem et al. [29] and Bert et al. [28]. The ST05 engines investigated in the first two experiments had slightly different dimensions and therefore the results corresponding to them should be treated rather as illustrative. The results from the work of Bert et al. [28] were obtained for the same engine unit as in the present simulations, except for the initial pressure load (pload =7bar), the measurements correspond to similar range of the rotational speeds and Th . The effect of the pload and rotational speed will be discussed later. From Fig. 11 it can be seen that the dependence of P on Th obtained in the present simulations and in [37, 29] is almost linear. According to [29] such a behaviour is typical and independent on the pload . However, in [29] it was found that the slope of the power decreases with the pload . Taking this into account the results obtained by Bert et al. [28] seem to be a bit surprising. It can be seen that in their

ACCEPTED MANUSCRIPT

case the indicated power diverges from the linear growth and moreover for the highest temperatures it is almost constant and equal to P ≈ 450W. These discrepancies could have been caused by measurement error, however, they should be kept in mind. Influence of the filling pressure on the indicated power is shown in Fig. 12. Consistently with [28, 29] the dependence P (pload ) is almost linear, the power increases with the filling pressure. Such a behaviour is rather expected since for the higher initial pressure the mass of the fluid inside the engine is larger and this ensures more effective heat transfer between the working spaces at the later stage. Unfortunately, in practice the initial filling pressure cannot be increased above a certain limit. It is bounded by the engine sealing and its strength potential to hold the gas at very high pressure. Comparing the results with the experimental data it can be seen that the agree ment in terms of the slopes is very good, particularly they overlaps with the data from work of Bert et al. [28]. However, it should be noted that their experiment had been performed for the higher Th value. The situation in which the experimental and numerical results fit well in some range of the filling pressures for different heater temperature is unrealistic. It can be attributed to the fact that in the simulations some losses have not been taken into account, i.e., a perfect mechanism and adiabatic walls without the heat losses assumption. Obviously, in experimental works such conditions cannot be achieved and some losses are unavoidable. This concerns apply to the results obtained by Hachem et al. [29], for exactly the same parameters as in our simulations (Th =773K, 600RPM). Compared to the present solutions the trend is exactly the same (linear) but the absolute values are smaller. The points marked in Fig. 12 are for the cases with the lower heater temperatures, Th = 673K and Th = 573K. They are shown here to illustrate that in our idealised conditions the obtained results cover the range of measured value for larger Th . Finally, we analyse dependence of the indicated power on the rotational speed. Figure 13 shows the plot of P (RPM) for Th =773K. It can be seen that an increase of the rotational speed leads to nearly linear increase of the power. In the experimental results such a behaviour is observed only for engine speeds smaller than 500RPM. The trends of present results agree relatively well with the measurements. Inconsistencies begin when the engine speeds up. In the works of Bert et al. [28] and

RI PT

Indicated Power [W]

600

400

SC

200

600

Exp. 7bar (Bert et al. (2014)) Exp. 10bar (Hachem et al. (2015)) Exp. 10bar (Alfarawi et al. (2016)) RANS 10bar 800

1000

M AN U

Heater Temperature [K]

AC C

EP

TE D

Figure 11: Indicated power as the function of heater temperature for 600RPM.

14

Hachem et al. [29] it was clearly shown that the engine has an optimal working speed for which it produces a maximum power. This optimal working point occurred at certain speed above which its further increase did not bring any profits, but rather lead to power decrease. The present numerical simulations were not able to capture such a point and moreover the results over-predicted the nominal power (Pmax = 500W) measured by the manufacturer (see Table 1). Most likely, this is because during the simulations the mechanical losses (e.g. friction losses) were not considered, and as it was noticed in [29], they are proportional to the average linear velocity of the pistons. Modelling of the mechanical losses in CFD is possible even for such complex object as the Stirling engine. This could be done by introducing additional terms into the flow equations, which role would be to mimic the pressure drop corresponding to the power needed to overcome the friction forces of the crank mechanism. Such an approach is planned for future studies and it is believed that it will improve the accuracy of the results. Summing up, the solutions presented in this subsection confirmed the validity of the applied solution algorithm. The results proved also that the CFD tools, even when used with some simplifications, can provide meaningful results.

ACCEPTED MANUSCRIPT

Exp. 973K (Bert et al. (2014)) Exp. 773K (Hachem et al. (2015)) RANS 773K

Exp. 7bar/773K (Bert et al. (2014)) Exp. 10bar/773K (Hachem et al. (2015)) RANS 10bar/773K

800

800

600

R600P10T673

400

600

RI PT

Indicated Power [W]

Indicated Power [W]

10bar/773K (Alfarawi et al. (2016))

400

R600P10T573

200

0

1

2

3

4

5

6

7

8

9

10

11

SC

200

12

200

600

800

1000

Engine speed [RPM]

Figure 12: Indicated power as the function of filling pressure for 600RPM.

4.3. Detailed flow characteristics

400

M AN U

Filling pressure [bar]

Figure 13: Indicated power over a range of engine speeds for Th =773K.

AC C

EP

TE D

In the current subsection we focus on detailed investigations of the flow dynamics inside the engine. This type of analysis would not be possible using simplified zero- or one-dimensional models. Here, we focus only on one configuration, i.e., R180 P10 T773 , which can be treated as an representative example for all the cases specified in Table 2. The aim of these investigations is to show the complexity of the flow phenomena in particular parts of the engine. We would like to emphasize the fact that the detailed flow analysis requires full 3D simulations without the geometrical simplifications (e.g. the assumption of cylindrical symmetry [33, 37, 48]) or simplifications related to the flow character (e.g. laminar regimes [33, 35, 24]), which may significantly affect the modelling accuracy. Sample flow patterns inside the cylinders, heater tubes and in the regions of the connecting tube inlets are demonstrated in Fig. 14. It shows isosurfaces of instantaneous vorticity magnitude at the time instances t0 , t0 + 10%Tcycle , t0 + 30%Tcycle and t0 + 50%Tcycle corresponding to different positions of displacer and power piston. The period of one full cycle is equal to Tcycle = 1/3s. Figure 14a corresponds to the time t0 , which has been chosen arbitrarily as a moment when the displacer moves towards the bottom dead point. According to the list of a-d characteristic positions discussed in Sec. 2.2, Fig. 14b corresponds to the situation

15

right after the point a), i.e., the time instance at t0 + 10%Tcycle , when the displacer starts moving up and the power piston is on the half-way down to the bottom dead location. The remaining figures cover later time instances of half of the cycle. The presented figures show that highest values of the vorticity occur in the heater tubes, particularly at the time instances when the displacer is far from the bottom dead point (Fig. 14a, c, d). This is related to the large values of the velocity of the air flowing between the regenerator and the cylinder above the displacer. On the contrary, when the displacer is close to the dead point the air velocity and its gradients are small, and thus the vorticity values in the tubes are significantly reduced, mainly in their upper parts. In general, in the heater tubes the smallest vorticity values occur in their curved part. However, the regions of the flow characterised by decreased vorticity do not occur in all tubes together at the same time but their locations change depending on the displacer position. It is somewhat surprising that the flow at the entrance to every single tube from/to the regenerator or from/to the cylinder has approximately the same features (velocity, turbulence intensity, etc.), which should lead to practically the same velocity/vorticity distributions inside the tubes. The analysis of the velocity field in the inlet tube regions has shown that such presumption is not fully correct. In fact, the flow veloc-

RI PT

ACCEPTED MANUSCRIPT

(b) t0 + 10%Tcycle

(d) t0 + 50%Tcycle

TE D

(c) t0 + 30%Tcycle

M AN U

SC

(a) t0

Figure 14: Isosurfaces of instantaneous vorticity magnitude equal to 500s−1 (blue), 1000s−1 (green), 1900s−1 (red).

AC C

EP

ity of the air entering the particular tube is largely conditioned on what happens in the cylinder. The presented figures show that the flow behaviour and vortical structures in this part of the engine are very complex, especially when the displacer moves down and the gas from the heater tubes is sucked to the cylinder (see Fig. 14a).

16

Besides the heater tubes the sub-figures in Fig. 14 show also the regions inside the cylinders and close to the endings of the tube connecting the cylinders. Evidently, the flow in these spaces seems to be much more complex than in the tubes. In general, the vorticity values are smaller but the shapes of its isosurfaces clearly indicate that the flow in these regions is non-uniform, non-symmetric and strongly turbulent. Figure 15 presents instantaneous temperature distributions in the symmetry plane and on the walls. The sub-figures show the results at the time instances denoted as t0 , t0 + 10%Tcycle , t0 + 20%Tcycle ,...., t0 + 50%Tcycle . The arrows indicate the direction in which the power piston and displacer move. It can be seen that when the displacer starts moving up (Fig. 15b) the hot air being in the cylinder space over the displacer, flows through the heater tubes and regenerator to the cooler surrounding the cylinder’s wall below the displacer. From the cooler, the air flows partially

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

(b) t0 + 10%Tcycle

(c) t0 + 20%Tcycle

AC C

EP

TE D

(a) t0

(d) t0 + 30%Tcycle

(e) t0 + 40%Tcycle

(f) t0 + 50%Tcycle

Figure 15: Instantaneous temperature distribution at the walls and the symmetry plane.

17

(a) t0

SC

RI PT

ACCEPTED MANUSCRIPT

(b) t0 + 10%Tcycle

(c) t0 + 30%Tcycle

(d) t0 + 50%Tcycle

(b) t0 + 10%Tcycle

AC C

EP

(a) t0

TE D

M AN U

Figure 16: Instantaneous velocity distribution (z-component) at the symmetry plane

(c) t0 + 30%Tcycle

(d) t0 + 50%Tcycle

Figure 17: Instantaneous streamlines of the velocity in the ”hot cylinder” and heater tubes coloured by z-velocity component.

18

to the tube connecting both cylinders and also to the space under the displacer. The fact that the flow in the heater changes its direction between the time moments t0 and t0 + 10%Tcycle is confirmed in Figs. 16a, b which show the contours of the vertical velocity in the main cross-section plane. The colours of the velocity contours in the heater tubes clearly indicate that the velocity changes the sign. It can be seen in Figs. 16b, c that in the regions connecting the tube with the cylinders the velocity is negative and quite large, approx. -7.5m/s. Closer inspection of the contours in the ”cold cylinder” shows that at its entrance the velocity is negative while on the opposite side (i.e. on the right side), close to the surface of the piston, it takes positive values. This reveals that when the power piston moves towards the bottom dead point the flow in the central part of the ”cold cylinder” recirculates counter-clockwise. Figures 15b,c,d show that the temperature distribution in the cooler space and inside the ”hot cylinder” under the displacer is non-symmetric. Significant differences in the air temperature in the heater tubes and in the region close to an entrance to the cooler show that the regenerator accumulates large amount of heat. Such a behaviour is desired and expected as the regenerator role is to store energy from the hot gas and to return it when the gas flows in the opposite direction. Figure 15e presents the results at the time instance when the power piston is in the bottom dead point and the displacer is in the middle location - this corresponds to the

ACCEPTED MANUSCRIPT

radially towards the centre and towards the outer walls. In effect two recirculation zones are created, i.e., large inner zone in the centre of the cylinder and the smaller one between the wall and the ”jets” of the gas flowing from the heater tubes. The occurrence of the recirculation zones is not obvious but it can be deduced from Fig. 16a, where the regions of positive and negative velocity are readily apparent. The fluid flow in the heater tubes and in the ”hot cylinder” in the space above the displacer is nicely presented in Fig. 17 showing the streamlines of instantaneous velocity coloured by the z-velocity component. Arrangement of the streamlines close to the surface of the displacer shows that the part of the air flows radially to the centre and also to the walls (see Fig. 17a). Some of the streamlines in the central part of the cylinder, and more precisely in the inner recirculation zone, are reminiscent of large distorted rings. When the displacer is close to the bottom dead point and its velocity is very small (Fig. 17b) the velocity of air is also small (green colour) and the flow behaves in a chaotic manner, which is reflected by irregular shapes of the streamlines. This situation changes when the displacer starts to move up quicker (Fig. 17c, d) and the compression process becomes stronger. It damps small scale motion and makes the flow uniform (see Fig. 16c, d). Lack of the velocity fluctuations above the displacer is manifested by the smooth and almost straight streamlines extending vertically across the whole cylinder space. The compression, forces the gas to flow through the heater tubes which total cross-section area is much smaller compared to the surface area of the displacer. As a result the flow inside the tubes strongly accelerates (red and blue colour).

EP

TE D

M AN U

SC

RI PT

characteristic point b) discussed in the Section 2.2. At this moment the cooler surrounding the lower part of the ”hot cylinder” is almost completely filled with the heated air. It can be seen that the hot gas is sucked to the space under the displacer through a small orifice between the walls connecting the cooler tube with casings of the ”hot cylinder” and cooler. When the power piston is in the bottom dead point the absolute values of the velocity in the ”cold cylinder” are close to zero (not shown here). The situation changes significantly when the piston starts moving back towards the top dead point. In this case the air being in the ”cold cylinder” flows to the cooler tube and forces the movement of the cold air to the ”hot cylinder”, i.e., in the space under the displacer. This can be directly observed in Fig. 16d but also it can be deduced from Fig. 15f. In the former case, it can be seen that the air entering the cylinder is oriented up (positive z-velocity), it flows through the cylinder space, hits on the bottom wall of the displacer and then moves along the wall and turns down near the cooler wall. Such a behaviour is a manifestation of existence the recirculation zone with the counter-clockwise rotation. Figure 15f shows that the cold air coming from the tube separates a stream of hot air flowing under the displacer from the cooler space. Although, the air from the tube flows mainly to the cylinder it also pushes on the gas being in the cooler. However, this kind of blockage occurs only locally in the place where the tube connects to the cooler. On the remaining part of its circumference the hot air gets to the cylinder through a narrow gap, as can be seen in Fig. 15f. Then it flows along the bottom wall towards the inner cylindrical wall and turns up. In effect, the flow in the ”hot cylinder” in the space under the displacer is not symmetric and exhibits strongly unsteady and turbulent behaviour. Compared to the previously described complex characteristic of the flow in the ”cold cylinder” the situation in the ”hot cylinder” in the region above the displacer seems to be much simpler. Figure 15 shows that at all time instances the temperature in this part is uniform and decreases slightly only when the displacer moves towards the top dead point. In the same time the velocity in the cylinder is small and almost uniform, as can be seen in Figs. 16c, d. The velocity field exhibits more complicated pattern only when the displacer is close to the bottom dead point or when it moves down. In this case the gas which is sucked from the heater tubes attaches to the displacer surface and flows

AC C

4.3.1. Characteristics of turbulent flow field Further differences in the flow behaviour when the displacer and the power piston bypass their bottom dead points are illustrated in Fig. 18. It presents the contours of the vorticity module (|Ω|), turbulent kinetic energy (TKE) and its dissipation rate () in the main cross-section plane at y = 0. Note that the last two quantities count in the overall budget of lost energy, they are related to inner flow losses. Unlike the mechanical losses related to an imperfection of the working elements (e.g. crank mechanism) or longitudinal and local viscous losses related to a pressure drop [29, 37], the nature of energy lost due to turbulence is different. The TKE 19

(b) |Ω| at t0 + 10%Tcycle

(e) TKE at t0

(f) TKE at t0 + 10%Tcycle

(c) |Ω| at t0 + 30%Tcycle

(d) |Ω| at t0 + 50%Tcycle

(g) TKE at t0 + 30%Tcycle

(h) TKE at t0 + 50%Tcycle

(k)  at t0 + 30%Tcycle

(l)  at t0 + 30%Tcycle

AC C

EP

TE D

M AN U

(a) |Ω| at t0

SC

RI PT

ACCEPTED MANUSCRIPT

(i)  at t0

(j)  at t0 + 10%Tcycle

Figure 18: Vorticity module (|Ω|), turbulent kinetic energy (TKE) and dissipation rate () showed at the symmetry plane.

20

ACCEPTED MANUSCRIPT

required to drive the flow at the assumed flow rate (Q), Ploss = Q × ∆p. The pressure distribution along the walls of the engine is shown in Fig. 19. Compared to ∆p ≈ 15kPa in the regenerator the pressure variations in the remaining parts are negligible. Turbulence is reborn only inside the tube connecting the cylinders. In addition, an increased level of TKE and  is visible under the displacer and inside the ”cold cylinder”. The occurrence of their largest values is directly related with the regions of large velocity (see Fig. 16). The maximum of TKE and  can be found in the upper part of the tube, where the flow issuing from the cooler hits the tube walls, and in the region right after the tube is sharply curved towards the ”cold cylinder”. Hence, one may assume that besides the heater tubes these are the places of the most intense production and dissipation of the turbulent kinetic energy, and hence, they significantly contribute to the drag forces. Finally, we analyse variations of the skin friction coefficient (f ) that directly affect the local and longitudinal pressure losses ∆pfrict ∝ f, µ, u. Its distributions along the engine walls in particular working phase are shown in Fig. 20. For example, at the time instances t0 and t0 + 10%Tcycle (Fig. 20a, b) the high values of f are located in parts next to the inlets of the tube connecting the cylinders and inside the ”cold cylinder”. The working fluid accelerates in those parts due to the suction caused by the motion of the power piston. Inversely, at the time instances t0 + 30%Tcycle and t0 + 50%Tcycle (Fig. 20c, d), the rapidly moving displacer imposes the acceleration of the working fluid inside the heater tubes. In this case the high values of f appear on the heater tube walls. Generally, it can be seen that the high values of f occupy regions where the working gas is accelerating as discussed previously (see Fig. 16). Acceleration of the flow is caused by the drastic pressure drop due to the substantial differences in the cross-sectional areas of the flow passages. The friction coefficient is directly proportional to the shear stress on the walls. Its growth is caused by the augmentation of the viscous stresses and it occurs in the regions of laminar-turbulent transition. The high values of f appearing on walls, e.g. at the connection tube at t0 + 10%Tcycle and heater tubes at t0 + 30%Tcycle , indicate that these parts are subjected to the high viscous losses. Additionally, the high values of  (see Figs. 18i-l) accompany the high values of f .

AC C

EP

TE D

M AN U

SC

RI PT

is contained in small scales and its production requires energy, which can be drawn only from the mean flow. In turn, the dissipation is responsible for transferring of TKE into heat through the viscous forces. A precise analysis of these two quantities and identification of the places where they can be reduced is unavoidable when using 1D models, it is possible only by CFD. Moreover, taking into account that the flow in the Stirling engine is strongly irregular and asymmetric, as shown above, the CFD should be performed for full 3D configurations without simplifying assumption, for instance, related to an azimuthal symmetry of the flow in the cylinders. The results presented in Fig. 18 show that in the ”hot cylinder” the flow exhibits turbulent behaviour only when the displacer moves down. In this case the values of both the vorticity and quantities characterising the turbulence level are large. Inside the heater tubes the distributions of TKE and  are typical for the pipe flows (i.e. small values are on the axes and large near the walls). On the other hand, in the region downstream the point where the tubes are connected to the cylinder both TKE and  show behaviour occurring in the jet flows. They are small in the potential core and large in the shear layers developing starting from the edges of tubes. It can be seen that these shear layers spread and meet together at the distance of approximately 3-4 tube diameters inside the cylinder. Close to this point the potential core collapses, the flow becomes fully turbulent and occupies significant part of the cylinder. The turbulent fluctuations in the ”hot cylinder” are quickly damped when the displacer is near the bottom dead point or when it moves up. It can be observed that at this stage the vorticity is large along the whole length of the tubes, whereas TKE and  are large mainly at their inlet parts and in the vicinity of the connection with the regenerator. These places can be regarded as the sources of significant energy losses. Because of large resistance of the regenerator the turbulence inside this part vanishes (TKE,  → 0) at a very close distance behind the tubes. The flow inside the regenerator and cooler becomes laminar (TKE,  ≈ 0) and remains almost uniform at all times (Fig. 16). However, it does not mean that the energy losses in the regenerator are small. In fact, as shown in [29] here they are the largest among all engine parts and most of them (44% [29]) are caused by the pressure drop (∆p) that has direct impact on the power (Ploss )

21

SC

RI PT

ACCEPTED MANUSCRIPT

(b) t0 + 10%Tcycle

M AN U

(a) t0

Figure 19: Absolute pressure distribution on the walls.

5. Conclusions

TE D

The later characterizes the rate at which the TKE is dissipated into the thermal energy.

The paper presented one of the first, multidimensional CFD study of the complete γ-type Stirling engine ST05G including its all geometrical details. The URANS simulations with realizable k −  and k − ω SST turbulence models were performed. The domain was discretised using moving grids employing layering and remeshing techniques to model the pistons movement. The comparison between the numerical and experimental results in terms of the indicated power as a function of engine rotational speed, filling pressure and heater temperature showed relatively good agreement. The obtained results have shown the same trends as in the experiments, thus confirming the qualitative correctness of the applied CFD model. Sensitivity analysis of the results concerning the regenerator parameters and turbulence modelling approach revealed their negligible impact on both the instantaneous solutions as well as the indicated engine power. The observed over-predictions of the exper-

(c) t0 + 30%Tcycle

(d) t0 + 50%Tcycle

EP

Figure 20: Skin friction coefficient distributions over the outer walls in particular working phases.

AC C

imental results were attributed to simplifying modelling assumption including the isothermal (heater and cooler) or perfectly adiabatic (all remaining) walls and lack of the mechanical losses caused by the kinematic mechanism. Certainly, they should be taken into account in future studies. Assuming the qualitative correctness of the applied model the main aim of this work was to capture and exhibit the complexity of the flow inside the Stirling engine including its transient features and existence of recirculation regions, separation points and streamline curvature - the phenomena characteristic for turbulent flows. It has been shown that the flow inside the engine has a 3D charac22

ACCEPTED MANUSCRIPT

SC

M AN U

Acknowledgements

EP

TE D

This work was supported by statutory funds BS/PB-1-103-3010/11/P. The first two authors are grateful to the Polish National Science Centre for funding their fellowships (Grant no: 2015/17/B/ST8/03217). PL-Grid infrastructure was used to carry out the computations. References

[7] Mardani A, Ghomshi AF. Experimental and numerical study of MILD combustion for gas turbine applications. Energy 2015;99:136-51. [8] Khalil AEE, Gupta AK. Swirling distributed combustion for clean energy conversion in gas turbine applications. Appl Energy 2011;88:3685-93. [9] Khalil AEE, Gupta AK. Distributed swirl combustion for gas turbine application. Appl Energy 2011;88:4898907. [10] Khalil AEE, Arghode VK, Gupta AK, Lee SC. Low calorific value fuelled distributed combustion with swirl for gas turbine applications. Appl Energy 2012;98:6978. [11] Tyliszczak A, Boguslawski A, Nowak D. Numerical simulations of combustion process in a gas turbine with a single and multi-point fuel injection system. Appl Energy 2016;174:153-65. [12] Arghode VK, Gupta AK. Investigation of reverse flow distributed combustion for gas turbine application. Appl Energy 2011;88:1096-104. [13] Arghode VK, Gupta AK, Bryden KM. High intensity colorless distributed combustion for ultra low emissions and enhanced performance. Appl Energy 2012;92:82230. [14] Arghode VK, Gupta AK. Investigation of forward flow distributed combustion for gas turbine application. Appl Energy 2011;88:29-40. [15] Arghode VK, Gupta AK. Development of high intensity CDC combustor for gas turbine engines. Appl Energy 2011;88:963-73. [16] Toftegaard MB, Brix J, Jensen PA, Glarborg P, Jensen AD. Oxy-fuel combustion of solid fuels. Prog Energ Combust 2010;36:581-625. [17] Scheffknecht G, Al-Makhadmeh L, Schnell U, Maier J. Oxy-fuel coal combustionA review of the current stateof-the-art. Int J Greenh Gas Con 2011;5:16-35. [18] Tyliszczak A, Rosiak A. LES-CMC simulations of a turbulent hydrogen jet in oxy-combustion regimes. Int J Hydrogen Energy 2016;41:9705-17. [19] Walker G. Stirling Engines. Oxford University Press; 1980. [20] Organ AJ. Stirling Cycle Engines: Inner Workings and Design. Wiley; 2007. [21] Goswami DY, Kreith F. Energy Conversion. CRC Press; 2007. [22] Kongtragool B, Wongwises S. A review of solar-powered Stirling engines and low temperature differential Stirling engines. Renew Sust Energ Rev 2003;7:131-54. [23] Kongtragool B, Wongwises S. A four power-piston lowtemperature differential Stirling engine using simulated solar energy as a heat source. Sol Energy 2008;82:493500. [24] Chen WL. A study on the effects of geometric parameters in a low-temperature-differential gamma-type Stirling engine using CFD. Int J Heat Mass Tran 2017;107:1002-13. [25] Thombare DG, Verma SK. Technological development in the Stirling cycle engines. Renew Sust Energ Rev 2008;12:1-38. [26] Parlak N, Wagner A, Elsner M, Soyhan HS. Thermodynamic analysis of a gamma type Stirling engine in nonideal adiabatic conditions. Renew Energ 2009;34:26673. [27] Timoumi Y, Tlili I, Nasrallah SB. Design and performance optimization of GPU-3 Stirling engines. Energy

RI PT

ter and therefore 1D or 2D analyses commonly applied in designing process cannot fully describe its physics. For instance, the analysis of the velocity field close to the heater pipes inlets revealed that the flow differs from pipe to pipe quite significantly. Hence, a usually made assumption that the flow behaviour in all tubes is identical is not fully correct. In most of the engine parts and particularly inside the hot cylinder or below the displacer the flow turned out to be very complex and characterized by small and large scale vortical structures. The temperature distributions within the cooler and inside the hot cylinder under the displacer were nonsymmetric. Analysis of the vorticity iso-surfaces and distributions of the turbulent kinetic energy has clearly indicated that the flow inside of almost all engine parts was spatially strongly non-uniform and turbulent. The applied CFD model allowed to reveal the regions characterized by large pressure drop and large values of the skin friction coefficient. It is believed that their precise identification can help to improve the overall engine efficiency.

AC C

[1] Flohr P, Stuttaford P. Combustors in gas turbine systems. In: Jansohn P, editor. Modern Gas Turbine Systems: High Efficiency, Low Emission, Fuel Flexible Power Generation. Woodhead Publishing; 2013. [2] G¨ uthe F, Hellat J, Flohr P. The Reheat Concept: The Proven Pathway to Ultralow Emissions and High Efficiency and Flexibility. ASME Turbo Expo. 2007;GT2007-27846. [3] G¨ uthe F, de la Cruz GM, Burdet A. Flue Gas Recirculation in Gas Turbine: Investigation of Combustion Reactivity and NOX Emission. ASME Turbo Expo. 200;,GT2009-59221. [4] Arghode VK, Gupta AK. Effect of flow field for colorless distributed combustion (CDC) for gas turbine combustion. Appl Energy 2010;87:1631-40. [5] Cavaliere A, de Joannon M. Mild Combustion. Prog Energ Combust 2004;30:329-66. [6] Kruse S, Kerschgens B, Berger L, Varea E, Pitsch H. Experimental and numerical study of MILD combustion for gas turbine applications. Appl Energy 2015;148:45665.

23

ACCEPTED MANUSCRIPT

2008;33:1100-14. [28] Bert J, Chrenko D, Sophy T, Le Moyne L, Sirot F. Simulation, experimental validation and kinematic optimization of a Stirling engine using air and helium. Energy 2014;78:701-12. [29] Hachem H, Gheith R, Aloui F, Nasrallah SB. Numerical characterization of a γ-Stirling engine considering losses and interaction between functioning parameters. Energ Convers Manage 2015;96:532-43. [30] Chen CL, Ho CE, Yau HT. Performance Analysis and Optimization of a Solar Powered Stirling Engine with Heat Transfer Considerations. Energies 2012;9:3573-85. ¨ ¨ Cetinkaya S, Saridemir S, Kara F. Ar[31] Ozgoren Y O, tificial neural network based modeling of performance of a beta type Stirling engine. Proc IMechE Part E: J Process Mechanical Engineering 2012;227:166-77. [32] Hooshang M, Moghadam RA, AlizadehNia S. Dynamic response simulation and experiment for gamma-type Stirling engine. Renew Energ 2016;86:192-205. [33] Leon-Salazar JL, Chen WL. A computational fluid dynamics study on the heat transfer characteristics of the working cycle of a β-type Stirling engine. Energ Convers Manage 2014;88:177-88. [34] Chen WL, Wong KL, Chang YF. A computational fluid dynamics study on the heat transfer characteristics of the working cycle of a low-temperature-differential γtype Stirling engine. Int J Heat Mass Tran 2014;75:14555. [35] Chen WL, Yang JL, Salazar YC. A CFD parametric study on the performance of a low-temperaturedifferential γ-type Stirling engine. Energ Convers Manage 2015;106:635-43. [36] Comsol multiphysics reference manual, https:// extras.csc.fi/math/comsol/3.4/doc/multiphysics/ wwhelp/wwhimpl/js/html/wwhelp.htm?context= multiphysics&file=html_commandintro.56.1.html; [accessed 13 August 2017]. [37] Alfarawi S, Al-Dadah R, Mahmoud S. Influence of phase angle and dead volume on gamma-type Stirling engine power using CFD simulation. Energ Convers Manage 2016;124:130-40. [38] Mahkamov K, Djumanov D. Three-dimensional CFD modeling of a Stirling engine. Proceedings of the 11th International Stirling Engine Conference 2003;97-107. [39] Mahkamov K. Design Improvements to a Biomass Stirling Engine Using Mathematical Analysis and 3D CFD Modeling. J Energy Resour Technol 2005;128:203-15. [40] Almajri AK, Mahmoud S, and Al-Dadah R. Modelling and parametric study of an efficient Alpha type Stirling engine performance based on 3D CFD analysis. Energ Convers Manage 2017;145:93-106. [41] Cheng CH, Chen YF. Numerical simulation of thermal and flow fields inside a 1-kW beta-type Stirling engine. Appl Therm Eng 2017;121:554-561. [42] The st05g stirling engine project. http: //ve-ingenieure.de/projekt_st05g_cnc_engl.html; [accessed 13 December 2016]. [43] Launder B, Spalding DB. Lectures in Mathematical Models of Turbulence. London: Academic Press; 1972. [44] Ansys Fluent 12.0 User’s Guide, http://users.ugent. be/~mvbelleg/flug-12-0.pdf; [accessed 20 July 2017]. [45] Pope SB. Turbulent Flows. Cambridge University Press; 2000. [46] Menter FR. Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications. AIAA J

AC C

EP

TE D

M AN U

SC

RI PT

1994;32:1598-1605. [47] Skaria M, Abdul Rasheed KK, Shafi KA, Kasthurirengan S, Upendra B. Computational Investigation on the performance of thermo-acoustically driven pulse tube refrigerator. IOP Conference Series: Materials Science and Engineering 2017;171:012087. [48] Mahkamov K. An Axisymmetric Computational Fluid Dynamics Approach to the Analysis of the Working Process of a Solar Stirling Engine. J Sol Energy Eng 2005;128:45-53. [49] Gheith R, Aloui F, Nasrallah SB. Study of the regenerator constituting material influence on a gamma type Stirling engine. J Mech Sci Technol 2012;26:1251-1255.

24

ACCEPTED MANUSCRIPT

The paper presents a multidimensional CFD study of the Gamma Stirling engine The simulations are performed using uRANS approach with k−epsilon and k−omega models The flow inside the engine has a 3D complex non-uniform character

AC C

EP

TE D

M AN U

SC

RI PT

The simulations reveal the regions characterized by strong vorticity values